% Test FFT with wavelet and added random noise
% TestFFTwithNoise.m  L Braile 3/18/06
npts = 200;   % Number of points in wavelet
dt = 0.01;    % sample interval (fnyq = 50 Hz)
freq = 5;     % set ~peak frequency of wavelet
timesh = 0.0; % allow timeshift to center the wavelet
nsw = 20;     % apply cosine bell taper to ends of signal
s = ricker(npts,freq,dt,timesh,nsw); % calculate ricker wavelet
t = [0:dt:(npts*dt - dt)];
noise = 0.1*randn(1,200); % calculate random noise
s = s + noise; % add random noise to data
figure
plot(t,s,'-r','linewidth',1.5)
set(gca,'fontsize',16,'linewidth',2)
xlabel('Time (s)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Ricker Wavelet','fontsize',16)
nf = 1024; % set length of fft
S = fft(s,nf); % calculate FFT of wavelet, the added zeros (nf = 1024)
% do not change the spectrum (except by providing finer sampling, df)
% because the zeros don't contribute to the sum of the area in the 
% Fourier integral
fnyq = 1/(2*dt); % Nyquist frequency
df = fnyq/(nf/2); % calculate frequency sample interval
f = [0:df:fnyq]; % calculate frequency variable (will be (nf/2) + 1 long
SS = S.*conj(S)/nf;  
SS = sqrt(SS); % calculate amplitude spectrum
figure
plot(f,SS(1:(nf/2)+1),'-r','linewidth',1.5)
set(gca,'fontsize',16,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Amplitude Spectrum of Ricker Wavelet','fontsize',16)
