fouriertimetest.asv

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

% Simple fourier filter coding for simulated time-based signal, where
% x = time in sec, millisec, or microsec. Original signal (blue)
% and filtered signal (red) are displayed superimposed in the top
% window, power spectrum and filter spectrum in bottom window.
% 'samplingtime' is the total duration of sampled signal, in seconds
%  (or millisec or microsec). 'samplerate' is the sample rate in Hz.
% 'centerfrequency' is the center frequency of filter, in Hz (or KHz or MHz).
% 'frequencywidth' is the frequency width of filter, in Hz (or KHz or MHz).
% 'Center' and 'width' are the center harmonic and width of
% the filter band in harmonics. '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.5, May, 2007

% Generate simulated signal
samplingtime=3; % Total duration of sampled signal, in sec, millisec, or microsec.
samplerate=100;  % Sample rate in Hz, KHz, or MHz.
n=samplingtime*samplerate;
frequency=20;
x=[0:(1/samplerate):samplingtime];
% Select one of these signal definitions by un-commenting it:
% y=sin(2*pi*frequency.*x); % Simple sine wave.
y=(1+sin(frequency./10*2*pi.*x)).*sin(frequency*2*pi.*x);  % f/10 Amplitude modulated sine wave
% y=gaussian(x,samplingtime/2,samplingtime/10); % Single gaussian band
% y=rand(size(x));  % Random white noise
% y=gaussian(x,samplingtime/2,samplingtime/10)+rand(size(x)); % Single gaussian band with random noise

% Initial values of filter parameters
centerfrequency=30;  % Center frequency of filter, in Hz.
center=centerfrequency*samplingtime; %  center harmonic (fourier component)
frequencywidth=1;  % Frequency width of filter, in Hz.
width=frequencywidth*samplingtime; %  width of filter (in harmonics)
shape=20; % 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)+1];
lft2=[(length(fy)/2):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(1:length(fy));  % 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    x=sec.')
subplot(2,1,2)
py=fy .* conj(fy); % Compute power spectrum
plotrange=1:length(fy)/2;
f=((plotrange-1)./samplingtime);
plot(f,real(py(plotrange)),f,max(real(py)).*ffilter(plotrange),'r')
title('BLUE = Power spectrum of signal      RED = Filter spectrum    x=Hz')
xlabel([ modestring '   Center = ' num2str(center) '    Width = ' num2str(width) '     Shape =  ' num2str(shape)])