Signal Processing

Adding Sines

% square wave sum of sins example

% change nfreqs to change accuracy of sq-wave approx.

n_long_periods = 2;
min_points_per_period = 128;
nfreqs = 256;
points_per_long_period = nfreqs * min_points_per_period;
npoints = n_long_periods * points_per_long_period;
points = linspace(1,npoints,npoints);
sum = zeros(1,npoints);
clf
figure
% Use only the odd frequencies
for freq = 1:2:nfreqs
newsin = sin(points*((2*pi*freq)/points_per_long_period));
wtsin = newsin / freq;
sum = sum + wtsin;

plot(wtsin, 'r');
hold on
end
plot(sum, 'b');
hold off

FFT Basics and watching a phasor

N = 64;
phase = 3;
periods = 4;
times= linspace(0,((N-1)/N)*2*pi*periods,N);
ptimes = times + (phase*2*pi*periods)/N;
anginc = (2*pi*4)/64;
signal1 = zeros(1,N);
signal2 = zeros(1,N);
for i = 1:N;
    signal1(i)=sin(times(i));
    signal2(i)=sin(ptimes(i));
end
plot(times,signal1, '-r');
hold on
plot(times,signal2, '-b');
y1 = fft(signal1,N);
y2 = fft(signal2,N);
pow1 = y1 .* conj(y1);
pow2 = y2 .* conj(y2);
figure
plot(1:N,pow1);
figure
plot(1:N,pow2);

Square wave as sum of sines example

% square wave sum of sins example

n_long_periods = 2;
min_points_per_period = 128;
nfreqs = 6;
points_per_long_period = nfreqs * min_points_per_period;
npoints = n_long_periods * points_per_long_period;
points = linspace(1,npoints,npoints);
sum = zeros(1,npoints);
clf
figure
% Use only the odd frequencies
for freq = 1:2:nfreqs
newsin = sin(points*((2*pi*freq)/points_per_long_period));
wtsin = newsin / freq;
sum = sum + wtsin;

plot(wtsin, 'r');
hold on
end
plot(sum, 'b');
hold off

Chirp (simple)...maybe look at autocorrelation?

y = zeros(1,100);
x = linspace(0, 6,100);
for i = 1:100
    y(i) = sin(x(i)^2.2);
end;
plot(x,y);

1-D deconvolution (with lots of debugging statements)

function ImOut = deconv1D(Im, mask,n);
% Blur and unblur a 1-d image w/ 1-D mask.
% here just imbedding inputs in 1xn vector, rest 0s.


plot(Im);
ImOut = zeros(1,n);

blurredim = conv (Im, mask)  % create blurred image

fftmask = fft(mask,n);
ifftmask = zeros(1,n);  % for 1/fft

fftIm = fft(blurredim,n);

for i = 1:n
ifftmask(i) = 1/fftmask(i);  % may create infsa
end

'ifftmask'
ifftmask
fixedifftmask = inffix(ifftmask,n); % remove infs
%fftmask .* ifftmask;
'should be ones'
fftmask .* fixedifftmask

'fftrow'
fftrow = fftIm .* fixedifftmask

%fftrow = inffix(fftrow,n);
ImOut = ifft(fftrow,n);
figure
plot(ImOut);
end

function fixout1 (X,n)
for i = 1:n
    if X(i) < 0
        X(i) = 0;
    end
end
end

function fixed = inffix(X,n)
fixed = X;
for i = 1:n
    if isinf(X(i))
        fixed(i) = 0.0;
 
    end
end
end

2-D deconvolution with 1-D blur function

function ImOut = deconvolve(Im, mask,n);
% Im an nxn image, mask the mask we want to de-blur by, n is
% a power of two.
%This is a cheat  since here I'm just doing one row at
% a time, since I know my mask is 1-D.

ImOut = zeros(n);

fftmask= fft(mask,n);
ifftmask = zeros(1,n);
for i = 1:n
ifftmask(i) = 1/fftmask(i);
end

ifftmask;
fixedifftmask = inffix(ifftmask,n);
%fftmask .* ifftmask;
'should be ones'
fftmask .* fixedifftmask;
size(fftmask);
size(ifftmask);

for i = 1:n
    fftrow = fft(Im(i,:));
size(fftrow);
fftrow = fftrow .* fixedifftmask;
fftrow = inffix(fftrow,n);
ImOutRow = ifft(fftrow,n);
ImOut(i,:) = ImOutRow;
end
ImOut = scale_it(ImOut,5.0);
end


function fixout1 (X,n)
for i = 1:n
    if X(i) < 0
        X(i) = 0;
    end
end
end

function fixed = inffix(X,n)
fixed = X;
for i = 1:n
    if isinf(X(i))
        fixed(i) = 0.0;
        i
    end
end
end

1-D Power Spectrum

function thePS = PowSpec1D(X,n)
 Y = fft(X,n);
 thePS = (Y .* conj(Y)) / n;
end

Blurring image with 1-D streak (motion blur)

function BlurData = make_data(InImage,xfilter)
 BlurData = conv2(InImage,xfilter);
 BlurData=scale_it(BlurData,1.0);
 imshow(BlurData);
end

Make the jet engine data

% .1 second of data with 3 sine waves in it.

wavecount = 3;
datalen = 128;
maxtime = .1;
times = linspace(0,.1,datalen);
waves = zeros(wavecount,datalen);
freqs = [60, 150, 350];
cycles = freqs*maxtime;
amps = [1, .4, 1.25];
mean = 0;
std = 1;  % or variance??

for i = 1:wavecount
   
    waves(i,:) = sin(times .*(2*pi*cycles(i))/maxtime);
end

noise = random('norm', mean, std, 1,datalen);

signal = sum(waves);
data = noise + signal;
%plot(times, signal);
%plot(times, noise);
plot(times,data);

Y= fft(data,datalen);
Pyy = fftshift(Y.*conj(Y))/datalen;
maxfreq = datalen/(2*maxtime);
freqaxis = linspace(-maxfreq, maxfreq,datalen);
figure
plot(freqaxis,Pyy);

Color Image to 2-D Power Spectrum Display

function sPS = PSPic(fn,N)
given a color image we think is about right size...
colarr = imread(fn);
%imshow(colarr);
bwarr = rgb_to_bw(colarr);

bwarr = double(colarr);
scaledbw = scale_it(bwarr, 1.00);
%imshow(scaledbw);
ft = fft2(scaledbw, N,N);
PS = real(ft .* conj(ft));
PS(1,1) = 0.0; % kill off central peak
PS = PS + 1.0; % make minimum value = 1
%PS = PS .^ (0.3); % OK but unusual
sPS = log(PS);    % works better
sPS = scale_it(sPS, 1.0);
figure;
imshow(fftshift(sPS));
% could also go back to 0-255 uint8s
ssps = uint8(sPS*255);
%figure;
%imshow(ssps);
end

function scaledX= scale_it( X, maxval)
%scale everything in X by division so that new maximum  is maxval
Xmax = max(max(X));
scaledX = X * (maxval/Xmax);
end

function bwarr = rgb_to_bw(arr)
% convert rgb to bw using some magic constants 
arr = double(arr);
bwarr = .3* arr(:,:,1) + .59*arr(:,:,2) + .11*arr(:,:,3);
end