Matlab / Octave scripts from Coursera' Computational Methods for Data Analysis week 1

fft01.m

% clear everything

clear all; close all; clc;

%create a gaussian

L=20
n=128;
x2=linspace(-L/2,L/2,n+1);
x=x2(1:n);
u=exp(-x.^2);
plot(x, u);

%perform fft

ut = fft(u);
plot(x, ut);

% get the gaussian back

plot(abs(fftshift(ut)));

% the 2 * pi / L is to get to rescale because Matlab thinks everything is 2*pi
% the next part is frequency numbers
% everything is shifted so we have to account that
% so we have to fft this as well

k=(2*pi/L) * [0:n/2-1 -n/2:-1];
plot(fftshift(k), abs(fftshift(ut)));

u = sech(x);
ud = -sech(x).*tanh(x);
u2d= sech(x) - 2*sech(x).^3;

ut=fft(u);

% ifft is the inverse fourier transform

uds = ifft((i*k) .*ut);
u2ds = ifft((i*k) .^2.*ut);
ks = fftshift(k);
plot(x, ud, 'r', x, uds, 'mo');

% if L is too small then we don't get a periodic function
% we have Gibbs oscillations at the end

fft02.m

% clear everything

clear all; close all; clc;

%create a gaussian

L=30
n=512;
t2=linspace(-L/2,L/2,n+1);
t=t2(1:n);
k=(2*pi/L) * [0:n/2-1 -n/2:-1];
ks=fftshift(k);

u=sech(t);
ut=fft(u);

% time
subplot(2,1,1), plot(t,u);

% spectrum
subplot(2,1,2), plot(ks, abs(fftshift(ut)));
axis([-25 25 0 1])

% if it's broad in time, narrow in spectrum and vice-versa

noise=1;
utn = ut+noise*(randn(1,n)+i*randn(1,n));
un=ifft(utn);
subplot(2,1,1), plot(t,u,'k',t, abs(un),'m');
subplot(2,1,2), plot(ks, abs(fftshift(ut))/max(abs(fftshift(ut))), 'k', ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm');

break

% ten times the noise

noise=10;

% note the i below means it's a complex number
utn = ut+noise*(randn(1,n)+i*randn(1,n));
un=ifft(utn);
subplot(2,1,1), plot(t,u,'k',t, abs(un),'m');
subplot(2,1,2), plot(ks, abs(fftshift(ut))/max(abs(fftshift(ut))), 'k', ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm');

% twenty times the noise
% without the original signal - can we still see it?

noise=10;
utn = ut+noise*(randn(1,n)+i*randn(1,n));
un=ifft(utn);

% use a filter
filter=exp(-k.^2);
utnf = filter.*utn;
unf = ifft(utnf);

% use a threshold to determine if there is a signal there

subplot(2,1,1), plot(t, abs(unf), 'g', t, 0 * t + 0.5, 'k');
subplot(2,1,2), plot(ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm', ks, fftshift(filter), 'b', ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm')

break

% move filter
filter=exp(-(k+15).^2);
utnf = filter.*utn;
unf = ifft(utnf);

% use a threshold to determine if there is a signal there

subplot(2,1,1), plot(t, abs(unf), 'g', t, 0 * t + 0.5, 'k');
axis([-15 15 0 1])
subplot(2,1,2), plot(ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm', ks, fftshift(filter), 'b', ks, abs(fftshift(utn))/max(abs(fftshift(utn))), 'm')

fft03.m

% clear everything

clear all; close all; clc;

%create a gaussian

T=30
n=512;
t2=linspace(-T/2,T/2,n+1);
% break it up into 512 points
t=t2(1:n);
% scaling
k=(2*pi/T) * [0:n/2-1 -n/2:-1];
% fft shifts things
ks=fftshift(k);

u=sech(t);
ut=fft(u);

noise=20;
utn=ut+noise*(randn(1,n)+i*randn(1,n));
un=ifft(utn);

subplot(2,1,1), plot(t, u, 'r', t, abs(un),'k')
subplot(2,1,2), plot(ks, abs(fftshift(ut)), 'r',ks,abs(fftshift(utn)))

break

# sampling

noise=30;
ave=zeros(1,n);
realizations=30;
for j=1:realizations
    utn=ut+noise*(randn(1,n)+i*randn(1,n));
    ave=ave+utn;
end
ave=abs(fftshift(ave))/realizations;
subplot(1,1,1)
plot(ks, abs(fftshift(ut)), 'r', ks, ave,'k', ks, abs(fftshift(utn)),'c')

fft04.m

% can't average moving signal in time domain
% but can do it in frequency domain 

% clear everything

clear all; close all; clc;

%create a gaussian

T=60
n=512;
t2=linspace(-T/2,T/2,n+1);
% break it up into 512 points
t=t2(1:n);
% scaling
k=(2*pi/T) * [0:n/2-1 -n/2:-1];
% fft shifts things
ks=fftshift(k);

slice =[0:0.5:10];
%deliberate mistake here because we redefine T
[T,S]=meshgrid(t, slice);
[X,S]=meshgrid(k, slice);
U=sech(T-10*sin(S)).*exp(i*0*T);
subplot(2,1,1)
waterfall(T,S,U), colormap([0 0 0]), view(-15,70)

% no waterfall plot in octave - see http://octave.1599824.n4.nabble.com/Waterfall-Plot-in-Octave-td4652110.html
% try
% surf(T,S,U)
% or
% plot3(T,S,U)
% it's not available in 3.6
% but maybe 3.8 - see http://www.gnu.org/software/octave/ and http://blogs.bu.edu/mhirsch/2013/12/compiling-octave-3-8/

for j=1:length(slice)
    UT(j,:)=abs(fftshift(fft(U(j,:))));
end

subplot(2,1,2)
% waterfall(fftshift(X),S,UT)
plot3(fftshift(X),S,UT) % no waterfall in octave

break

% shift the frequency by 10

U=sech(T-10*sin(S)).*exp(i*10*T);
for j=1:length(slice)
    UT(j,:)=abs(fftshift(fft(U(j,:))));
end

subplot(2,1,1)
U=real(U);
plot3(T,S,abs(U))
UT=real(UT);    % required by octave - see http://ubuntuforums.org/showthread.php?t=1025604
subplot(2,1,2)
plot3(fftshift(X),S,UT)

break

% get rid of the phase

subplot(2,1,1)
U=real(U);
plot3(T,S,abs(U))

break

% add noise
noise=20;
for j=1:length(slice)
    UT(j,:)=abs(fftshift(fft(U(j,:))+noise*(randn(1,n)+i*randn(1,n))));
end
UT=real(UT);

subplot(2,1,1)
U=real(U);
% waterfall(T, S, abs(U))
plot3(T,S,abs(U)) % no waterfall in octave
UT=real(UT);    % required by octave
subplot(2,1,2)
plot3(fftshift(X),S,UT)

break

% what does the noise look like in the time domain?

noise=20;
for j=1:length(slice)
    UT(j,:)=fft(U(j,:))+noise*(randn(1,n)+i*randn(1,n));
    UN(j,:)=ifft(UT(j,:));
end
UT=real(UT);

subplot(2,1,1)
plot3(T,S,abs(real(UN)))

This is a rough transcription - feedback, things I have missed etc are welcome.