% TestMovingAverageFilter.m
% L Braile 4/6/05
% simple wavelet
temp = [0 0.3 1.5 3 0 -3 -2 0.8 0.8 0.1 0]; % coarse sample interval
nt = length(temp);
n4 = 2*nt;
dt = 0.5;
t1 = [0:1:nt-1];
t2 = [0:dt:nt-1];
s4 = spline(t1,temp,t2); % interpolate to ~ twice as many samples
% simple wavelet "buried" in noise (begins at location 50)
s5 = randn(1,100); % make random noise scaled to std of 1
nt = length(s4);
s5(50:50+nt-1) = s5(50:50+nt-1) + s4(1:nt); % add signal to noise
n5 = length(s5);
nf = 5; % length of filter, odd number
f = ones(1,nf)/nf; % design moving average filter
s = conv(s5,f); % Apply moving average filter with convolution
nh = (nf-1)/2;
t = [0:1:n5-1];
sf = s(1,nh+1:n5+nh); % delete points at beginning and end
% of filtered output to align with input signal
plot(t,s5,'-b') % Plot unfiltered signal
hold on
plot(t,sf,'-r','linewidth',3) % Plot filtered signal
set(gca,'fontsize',16,'linewidth',2)
xlabel('Sample Number','fontsize',16)
ylabel('Amplitude','fontsize',16)
title('Moving Average Filter: Thin line = unfiltered, Bold line = filtered','fontsize',16)
hold off
% calculate and plot FFT of signal before and after filtering
fnyq = 1/(2*dt);
nfreq = 1024;
df = fnyq/(nfreq/2);
freq = [0:df:fnyq]; % define nfreq/2 + 1 length array for plotting
                    % spectra (0 to Nyquist Frequency)
X1 = fft(s5,nfreq);  % unfiltered spectrum
X2 = fft(sf,nfreq);  % Filtered spectrum
x1 = sqrt(X1.*conj(X1)); % amplitude spectrum
x2 = sqrt(X2.*conj(X2));
figure
loglog(freq,x1(1:(nfreq/2)+1),'b--',...
    freq,x2(1:(nfreq/2)+1),'r-','linewidth',1.5)
set(gca,'fontsize',16,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('Relative Amplitude','fontsize',16)
title('Amplitude Spectrum of Unfiltered (--) and Filtered (-) Signals','fontsize',16)
axis([0.01 1 0.1 100]);