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
Last update: 04/22/2011: RN