%  Transform Tool  (L Braile, October 2001; revised January 2002)
% Calculate Fourier Transforms os selected signals and display  
% Amplitude and Phase Spectra in Linear, Log and Semi-Log axes
% define signals each 201 points long, dt=0.01
dt=0.01;
t=0:dt:2;
%
% b1 short boxcar
b1=[zeros(1,98) ones(1,5) zeros(1,98)];
% b2  boxcar
b2=[zeros(1,95) ones(1,11) zeros(1,95)];
% b3  boxcar
b3=[zeros(1,90) ones(1,21) zeros(1,90)];
% b4  boxcar
b4=[zeros(1,75) ones(1,51) zeros(1,75)];
% b5  boxcar
b5=[zeros(1,50) ones(1,101) zeros(1,50)];
% sq square wave
sq=[zeros(1,75) ones(1,25) -1*ones(1,25) zeros(1,76)];
% sp spike
sp=[zeros(1,100) 1 zeros(1,100)];
% st step function
st=[zeros(1,50) ones(1,151)];
% ra random number
randn('state',1);
ra=randn(1,201);
ramax=max(abs(ra));
ra=ra/ramax;
% un uniformly distributed random number
un=rand(1,201);
un=2*un - 1;
% c1 low frequency cosine
%c1=cos(2*pi*1*t);
% c2 high frequency cosine
c2=cos(2*pi*5*t);
% Sum of 3, 8, 20 Hz sins and noise
sn=sin(2*pi*3*t)+sin(2*pi*8*t)+sin(2*pi*20*t)+ra;
smax=max(abs(sn));
sn=sn/smax;
% si sin x/x or sinc function
sit=sin(2*pi*4*t(2:101))./(2*pi*4*t(2:101));
si=[fliplr(sit) 1 sit];
% df differentiation operator
df=[zeros(1,99) 1 -1 zeros(1,100)];
% t1 short triangle
tr1=[1:1:5]/5;
tr2=[4:-1:1]/5;
t1=[zeros(1,96) tr1 tr2 zeros(1,96)];
% t2 short triangle
tr1=[1:1:26]/26;
tr2=[25:-1:1]/26;
t2=[zeros(1,75) tr1 tr2 zeros(1,75)];
% t3 long triangle
tr1=[1:1:51]/51;
tr2=[50:-1:1]/51;
t3=[zeros(1,50) tr1 tr2 zeros(1,50)];
% ga, gaussian
ga1=exp(-(32*t(2:101).*t(2:101)));
ga=[fliplr(ga1) 1 ga1];
% cb cosine bell
nsw=100;
   for i = 1:nsw
      j = nsw-i;
      fa1(i) = 0.5*(cos(pi*j/nsw) + 1);
   end
cb=[fa1 1 fliplr(fa1)];
% ex exponential
thalf=0.25;
x1=t(2:151)*0.696/thalf;
ex1=exp(-x1);
ex=[zeros(1,50) 1 ex1];
% w1 Ricker wavelet
%[rk]=ricker(201,20,dt,0,25);
%rmax=max(abs(rk));
%w1=rk/rmax;
% w2 Ricker wavelet
[rk]=ricker(201,5,dt,0,25);
rmax=max(abs(rk));
w2=rk/rmax;
% s1 spike series 1
s1=[zeros(1,10) 1 zeros(1,80) -0.7 zeros(1,100) 0.4 zeros(1,8)];
% s2 spike series 2
s2=[zeros(1,10) 1 zeros(1,45) .4 zeros(1,34) -0.7 zeros(1,55) ...
      -.3 zeros(1,44) 0.4 zeros(1,8)];
% s3 spike series 3
s3=[zeros(1,50) 1 zeros(1,8) -0.2 zeros(1,9) 0.7 zeros(1,8) -0.8 zeros(1,5) 0.3...
      zeros(1,12) -0.4 zeros(1,8) 0.25 zeros(1,5) 0.1 zeros(1,9) -0.1 zeros(1,6)...
      0.1 zeros(1,71)];
% rp ramp
%rp1=[1:1:101]/101;
%rp=[zeros(1,50) rp1 zeros(1,50)];
% sw sweep 1-5 Hz, 1000 points, dt=0.01
for i=1:183
   ti=(i-1)*dt;
   f=2 + 1*(ti);
   sw1(i)=sin(2*pi*f*ti);
end
sw=[zeros(1,9) sw1 zeros(1,9)];
% di Dirac comb
di=[zeros(1,10) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 ...
   zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 zeros(1,19) 1 ...
   zeros(1,10)];
% w3 wavelet - signal n=2
ts=0.21; n=2; npts=201; na=25; thalf=0;
[w3]=signal(ts,dt,n,0,npts,na,thalf);
% w4 wavelet - signal n=4
%ts=0.41; n=4; npts=201; na=25; thalf=0;
%[w4]=signal(ts,dt,n,0,npts,na,thalf);
% w5 wavelet - signal n=5
ts=0.51; n=5; npts=201; na=25; thalf=0;
[w5]=signal(ts,dt,n,0,npts,na,thalf);
% w6 wavelet - signal n=7
ts=0.71; n=7; npts=201; na=25; thalf=0;
[w6]=signal(ts,dt,n,0,npts,na,thalf);
% w7 wavelet - signal n=7 thalf=0.15
%ts=0.71; n=7; npts=201; na=25; thalf=0.15;
%[w7]=signal(ts,dt,n,0,npts,na,thalf);
% w8 wavelet - signal n=9 thalf=0.2
ts=0.91; n=9; npts=201; na=25; thalf=0.2;
[w8]=signal(ts,dt,n,0,npts,na,thalf);
%
%
% plot signals
figure
subplot(5,6,1), plot(t,b1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'b1, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
%
subplot(5,6,2), plot(t,b2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'b2, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,3), plot(t,b3,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
title('                              Functions and Signals for Fourier Transform Tool',...
   'fontweight','bold','fontsize',12)
text(0.1,-0.9,'b3, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,4), plot(t,b4,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'b4, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,5), plot(t,b5,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'b5, boxcar','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,6), plot(t,sq,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sq, sq. wave','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,7), plot(t,sp,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sp, spike','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
%
subplot(5,6,8), plot(t,df,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'df, differentiator','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,9), plot(t,st,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'st, step func.','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,10), plot(t,ra,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ra, random nor.','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,11), plot(t,un,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'un, random unif.','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,12), plot(t,c2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'c2, cosine','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,13), plot(t,si,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'si, sinc','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
ylabel('Relative Amplitude','FontWeight','bold')
%
subplot(5,6,14), plot(t,ga,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ga, gaussian','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,15), plot(t,t1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'t1, triangle','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,16), plot(t,t2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'t2, triangle','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,17), plot(t,t3,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'t3 triangle','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,18), plot(t,cb,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'cb, cosine bell','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,19), plot(t,ex,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'ex, exponential','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
%
subplot(5,6,20), plot(t,di,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'di dirac comb','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,21), plot(t,s1,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'s1 spike series','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,22), plot(t,s2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'s2, spike series','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,23), plot(t,s3,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'s3, spike series','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,24), plot(t,sw,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sw, sweep','FontWeight','bold')
set(gca,'XTickLabelMode','manual','XTickLabel',[])
set(gca,'YTickLabelMode','manual','YTickLabel',[])
%
subplot(5,6,25), plot(t,sn,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'sn, 3 sines+ra','FontWeight','bold')
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,26), plot(t,w2,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w2, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,27), plot(t,w3,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w3, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,28), plot(t,w5,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w5, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,29), plot(t,w6,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w6, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
subplot(5,6,30), plot(t,w8,'r-','linewidth', 1.5), axis([0 2 -1.1 1.1])
text(0.1,-0.9,'w8, wavelet','FontWeight','bold')
set(gca,'YTickLabelMode','manual','YTickLabel',[])
xlabel('Time (s)','FontWeight','bold')
%
%
pause
%
%
while 1
%
sig1=input(...
   'Type 2 character name for signal to be Fourier Transformed or qq to quit ');
%
if strcmp(sig1,'qq') | strcmp(sig1,'QQ'), break, end
%
% Plot signal
%
figure
subplot(2,2,1), plot(t,sig1,'r-','linewidth',2), axis([0 2 -1.1 1.1])
title('Signal','FontWeight','bold')
xlabel('Time (s)','FontWeight','bold')
ylabel('Relative Amplitude','FontWeight','bold')
%
% Calculate Fourier Transform and Amplitude and Phase Spectra
n=1024;
s=sig1;
S=fft(s,n);
fnyq=1/(2*dt);
freq = (0:n/2)/(n/2)*fnyq;
% Apply phase shift to time series (in Freq domain) to make 0 - 2 sec signals centered
% on origin -- -1 to +1 seconds.
%tshift=1;  % shift is 1 second
%theta=2j*pi*freq*tshift;
%e=exp(theta);
%S(1:n/2+1)=S(1:n/2+1).*e;
%End Phase shift adjustment
M=abs(S);
P=unwrap(angle(S));
%P=angle(S);
%
% Linear - Linear Plot
subplot(2,2,2)
[haxes,hline1,hline2] = plotyy(freq,M(1:(n/2)+1),...
   freq,P(1:(n/2)+1)*180/pi);
title('Amplitude and Phase Spectra -- Linear Axes','FontWeight','bold','fontsize',10)
set(haxes,'fontsize',10)
set(hline1,'color','r','linewidth',2)
set(hline2,'color','b','linewidth',2)
xlabel('Frequency (Hz)','fontsize',10,'FontWeight','bold')
axes(haxes(1))
v=axis;
axis([0 50 v(3) v(4)]);
ylabel('Amplitude Spectrum (Relative Ampl.)','color','r','FontWeight','bold','fontsize',10)
axes(haxes(2))
v=axis;
axis([0 50 v(3) v(4)]);
ylabel('Phase Spectrum (Degrees)','color','b','FontWeight','bold','fontsize',10)
%
% Log - Log Plot
subplot(2,2,3)
[haxes,hline1,hline2] = plotyy(freq,M(1:(n/2)+1),...
   freq,P(1:(n/2)+1)*180/pi,'loglog','loglog');
title('Amplitude and Phase Spectra -- Log Axes','FontWeight','bold','fontsize',10)
%set(haxes,'fontsize',10)
set(hline1,'color','r','linewidth',2)
set(hline2,'color','b','linewidth',2)
axes(haxes(1))
v=axis;
axis([0.1 100 v(3) v(4)]);
hold on
set(gca,'XAxisLocation','top','XTickLabel',[])
hold off
ylabel('Amplitude Spectrum (Relative Ampl.)','color','r','fontsize',10,'FontWeight','bold')
axes(haxes(2))
xlabel('Frequency (Hz)','fontsize',10,'FontWeight','bold')
v=axis;
axis([0.1 100 v(3) v(4)]);
ylabel('Phase Spectrum (Degrees)','color','b','fontsize',10,'FontWeight','bold')
%
% Semi-Log Plot
subplot(2,2,4)
[haxes,hline1,hline2] = plotyy(freq,M(1:(n/2)+1),...
   freq,P(1:(n/2)+1)*180/pi,'semilogy','semilogy');
title('Amplitude and Phase Spectra -- Semi-Log Axes','FontWeight','bold','fontsize',10)
%set(haxes,'fontsize',10)
set(hline1,'color','r','linewidth',2)
set(hline2,'color','b','linewidth',2)
axes(haxes(1))
v=axis;
axis([0 50 v(3) v(4)]);
hold on
set(gca,'XAxisLocation','top','XTickLabel',[])
hold off
ylabel('Amplitude Spectrum (Relative Ampl.)','color','r','fontsize',10,'FontWeight','bold')
axes(haxes(2))
xlabel('Frequency (Hz)','fontsize',10,'FontWeight','bold')
v=axis;
axis([0 50 v(3) v(4)]);
ylabel('Phase Spectrum (Degrees)','color','b','fontsize',10,'FontWeight','bold')
%
end