1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
| % 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:dtnpts*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,'linewid th',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(1nf/2)+1),'-r','linewidth',1.5)
set(gca,'fontsize',16,'linewid th',2)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Amplitude','fontsize', 16)
title('Amplitude Spectrum of Ricker Wavelet','fontsize',16) |
