sampling

% This is our CT function
clear
close all
func = @(t) -gaus(1.5*t+1)+ 0.8*gaus(t-0.3) -0.6*gaus(3*t) + gaus(t+0.7)+0.6*gaus(1.5*t-1.5) - 0.3*gaus(t/3);

tn = linspace(-4,4,501);
yn = func(tn);
plot(tn,yn,'LineWidth',2);

% Real part of Fourier transform (even symmetry)
fourier(func,3,3);

maximum real value 0.486768
maximum imag value 0.489684

% get Fourier transform, show magnitude (absolute value)
close all;
N = 480;
t = (0:N-1)*24/N - 12;
dt = t(2)-t(1);
x = func(t);
x = fftshift(x);
[fx ft] = fourier_xform(x,t);
fa = abs(fx)*dt;
idx = find(abs(ft)<4);
plot(ft(idx),fa(idx),'LineWidth',2);
%plot(ft(idx),real(fx(idx)),ft(idx),imag(fx(idx)),'LineWidth',2);
xlabel('frequency');
ylabel('FT magnitude');

% calculate some basic parameters
fs= 1/dt;
fprintf('sample frequency %g\n',fs);
N = length(t);
period = N*dt;
fprintf('sample time %g\n',period);
f0 = 1/period;
fprintf('FFT fundamental %g\n',f0);
fmax = max(ft);
fprintf('maximum FFT frequency %g\n',fmax);

sample frequency 20
sample time 24
FFT fundamental 0.0416667
maximum FFT frequency 9.95833

% show sprectrum for sampled waveform with specified sampling frequency
fs = 8;
freq = linspace(-8,8,801);
y = zeros(size(freq));
for k=-5:5
    y = y + pfunc(freq,fa,ft+fs*k);
end
plot(freq,y,'LineWidth',2);
xlabel('frequency');
ylabel('FT magnitude');

% show samples for sampled waveform
fs = 8;
N = 8*fs + 1;
if N%2 == 0
    N=N+1;
end
ts = linspace(-4,4,N);
ys = pfunc(ts,yn,tn);
stem(ts,ys);
xlabel('time');
ylabel('signal level');
hold on;
plot(tn,yn,'r--');
hold off;

% generate a movie showing the effects of decreasing the sample frequency
aviobj = avifile('test1.avi','compression','Cinepak');
fig = figure;
for fs = linspace(12,2,120);

    y = zeros(size(freq));
    for k=-5:5
        y = y + pfunc(freq,fa,ft+fs*k);
    end
    plot(freq,y,'LineWidth',2);
    axis([-8 8 0 0.8]);
    xlabel('frequency');
    ylabel('FT magnitude');
    F = getframe(fig);
    aviobj = addframe(aviobj,F);
end
close(fig);
aviobj = close(aviobj);

% generate a movie showing the effects of decreasing the sample frequency
% version 2
aviobj = avifile('test2.avi','compression','Cinepak');
freq = linspace(-8,8,801);
fig = figure;
for fs = linspace(10,1,120);
    subplot(2,1,1);
    N = 81;
    ts = linspace(-4,4,N)*10/fs;
    ts = ts(abs(ts)<=4);
    ys = pfunc(ts,yn,tn);
    stem(ts,ys);
    xlabel('time');
    ylabel('signal level');
    hold on;
    plot(tn,yn,'r--');
    hold off;
    subplot(2,1,2);
    y = zeros(size(freq));
    for k=-5:5
        y = y + pfunc(freq,fa,ft+fs*k);
    end
    plot(freq,y,'LineWidth',2);
    axis([-8 8 0 1]);
    xlabel('frequency');
    ylabel('FT magnitude');
    F = getframe(fig);
    aviobj = addframe(aviobj,F);
end
close(fig);
aviobj = close(aviobj);