fouriertest.asv

A fourier filter for time-series signals. Does not require the Signal Processing Toolbox. There is

% Simple fourier filter coding for simulated time-series signal,
% without the interactive sliders. Original signal (blue)
% and filtered signal (red) are displayed superimposed in the top
%  window, power spectrum and filter spectrum in bottom window.
% 'Center' and 'width' are the center harmonic and width of
% the filter band. 'Shape' determines the sharpness of the 
% cut-off. If shape = 1, the filter is Gaussian; as shape 
% increases the filter shape becomes more and more 
% rectangular. Set mode = 0 for band-pass filter, 
% mode = 1 for band-reject (notch) filter.
%  T. C. O'Haver ([e-mail]), version 1.4, April, 2007

% Generate simulated signal
n=400;
x=[1:4:n];
% y=sin(2*pi.*x./n); % Simple sine wave with frequency = 1.
% y=(1+sin(10*pi.*x./n)).*sin(100*pi.*x./n);  % Amplitude modulated sine wave
y=sin(20*pi.*x./n)+sin(40*pi.*x./n)+sin(60*pi.*x./n)+5.*gaussian(x,length(x)/2,length(x)/10); % Gaussian abnd with three superimposed sine waves
% y=gaussian(x,n/2,n/10); % Single gaussian band
% y=rand(size(x));  % Random white noise

% Initial values of filter parameters
center=1; %  center harmonic (fourier component)
width=1; %  width of filter (in harmonics)
shape=2; % filter shape (higher numbers = sharper cutoff)
mode=0;  % mode=0 for band-pass filter, mode=1 for band-reject (notch) filter

% Fourier filter 
fy=fft(y); % Compute Fourier transform of signal y
% Compute filter shape
lft1=[1:(length(fy)/2)];
lft2=[(length(fy)/2+1):length(fy)];
ffilter1=ngaussian(lft1,center+1,width,shape);
ffilter2=ngaussian(lft2,length(fy)-center+1,width,shape);
ffilter=[ffilter1,ffilter2];
modestring='Band-pass mode:';
if mode==1, 
    ffilter=1-ffilter; 
    modestring='Band-reject (notch) mode:';
end
if length(fy)>length(ffilter), ffilter=[ffilter ffilter(1)];,end
ffy=fy.*ffilter;  % Multiply filter by Fourier transform of signal
ry=real(ifft(ffy)); % Inverse transform to recover filtered signal 'ry'

% Plot original and filtered signal in top plot, 
% power spectrum and filter in lower plot
subplot(2,1,1)
plot(x,y,x,ry,'r')
title('BLUE = Original signal        RED = Filtered signal')
subplot(2,1,2)
py=fy .* conj(fy); % Compute power spectrum
plotrange=1:length(fy)/2;
plot(plotrange-1,real(py(plotrange)),plotrange-1,max(real(py)).*ffilter(plotrange),'r')
title('BLUE = Power spectrum of signal      RED = Filter spectrum')
xlabel([ modestring '   Center = ' num2str(center) '    Width = ' num2str(width) '     Shape =  ' num2str(shape)])