% eda06_12
% makes exemplary time series d(t),
% amplitude spectral density s(f) pairs for
% Crib Sheet 6.2

clear all;

% setup for time
N=128;
Dt = 0.5;
t = Dt*[0:N-1]';

% setup for Fourier transform
T=N*Dt;
fmax=1/(2*Dt);
Df =fmax/(N/2);
Nf=N/2+1;
f=Df*[0:N/2,-N/2+1:-1]';

figure(1);
clf;
sp=1;

% narrow function
sigmad = 1.0;
d = exp(-((t-T/2).^2)/(2*(sigmad^2))) + random('Normal',0,0.1,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, 0, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

% zero area wavelet
sigmad1 = 2.5;
sigmad2 = 1.0;
d = exp(-((t-T/2).^2)/(2*(sigmad1^2)));
d = d.*cos((t-T/2)/sigmad2);
d = d + random('Normal',0,0.1,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

% wide function
sigmad = 7.0;
d = exp(-((t-T/2).^2)/(2*(sigmad^2))) + random('Normal',0,0.1,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, 0, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

% triangle
d = [zeros(N/4,1); (0:N/4-1)'/(N/4); (N/4-1:-1:0)'/(N/4); zeros(N/4,1)];
d = d + random('Normal',0.0,0.005,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

figure(2);
clf;
sp=1;

% approximately sinusoidal
w0 = 0.25*2*pi*(0.95+0.1*[0:N-1]'/(N-1));
d = sin( w0 .* t ) + random('Normal',0,0.2,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

% beating
w1 = 0.25*2*pi;
w2 = w1/2;
d = sin( w1*t ) + 0.5*sin( w2*t ) + random('Normal',0,0.2,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

% chirp
w0 = 0.25*2*pi*(2*[0:N-1]'/(N-1));
d = sin( w0 .* t ) + random('Normal',0,0.05,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);

% sawtooth
d = mod(t,N/10)/(N/10)-0.5;
d = d + random('Normal',0.0,0.005,N,1);
dtilde = Dt * fft(d);
s2=(2/T)*abs(dtilde).^2;
s=sqrt(s2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('t');
ylabel('d(t)');
axis( [0, T, -1.1, 1.1] );
plot(t,d,'k-','LineWidth',2);

subplot( 4, 2, sp);
sp=sp+1;
set(gca,'LineWidth',2);
hold on;
xlabel('f');
ylabel('s(f)');
axis( [0, fmax, 0, max(s)] );
plot(f(1:Nf),s(1:Nf),'k-','LineWidth',2);





