redrawfourierfilter.m

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

function ry=RedrawFourierFilter(xvector,yvector,centerfrequency,filterwidth,filtershape,mode)
% Separate graph windows for the original and filtered signals.
% Computes and plots fourier filter for signal yvector.  Centerfrequency
% and filterwidth are the center frequency and width of the
% pass band, in harmonics. 'filtershape' determines the sharpness of the 
% cut-off. If filtershape =1, the filter is Gaussian; as filtershape 
% increases the filter shape becomes more and more rectangular. 
% mode=0 for band-pass filter, linear plot; mode=1 for  band-pass filter,
% semilogx plot; mode=2 for band-reject (notch) filter, linear plot;
% mode=3 for band-reject (notch) filter, semilogx plot.
% In this version, the x-axis of the bottom window is expressed in Hz
% and CenterF and WidthF are the filter center and width in Hz.
%  T. C. O'Haver ([e-mail]),  version 1.5, May, 2007
global signalstring
fy=fft(yvector);
lft1=[1:(length(fy)/2)];
lft2=[(length(fy)/2+1):length(fy)];
% Compute filter shape
ffilter1=ngaussian(lft1,centerfrequency+1,filterwidth,filtershape);
ffilter2=ngaussian(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
ffilter=[ffilter1,ffilter2];
modestring='Band-pass:';
if mode>1, 
    ffilter=1-ffilter; 
    modestring='Band-reject (notch):';
end
if length(fy)>length(ffilter), ffilter=[ffilter ffilter(1)];,end
ffy=fy.*ffilter;  % Multiply filter by Fourier Transfork of signal
ry=real(ifft(ffy));
% Plot original and filtered signal in top plot, power spectrum and filter
% in lower plot
subplot(3,1,1)
plot(xvector,yvector,'b');
title(['Original Signal = ' num2str(signalstring) '    x=sec.'])
axis([xvector(1) xvector(length(xvector)) min(yvector) max(yvector)]);  

subplot(3,1,2)
plot(xvector,ry,'r');
title('Filtered signal    x=sec.')
axis([xvector(1) xvector(length(xvector)) min(yvector) max(yvector)]);  

subplot(3,1,3)
py=fy .* conj(fy); % Compute power spectrum
plotrange=1:length(fy)/2;
f=((plotrange-1)./range(xvector));
switch mode,
  case {0,2}
    plot(f,real(py(plotrange)),f,max(real(py)).*ffilter(plotrange),'r')
  case {1,3}
    semilogx(f,real(py(plotrange)),f,max(real(py)).*ffilter(plotrange),'r')
  otherwise,
end
title('BLUE = Power spectrum of signal      RED = Filter spectrum   x=Hz')
xlabel([ modestring '   CenterFreq = ' num2str(centerfrequency./range(xvector)) '    FreqWidth = ' num2str(filterwidth./range(xvector)) '     Shape =  ' num2str(filtershape)])
if centerfrequency+filterwidth/2<length(fy)/4,
    axis([0 max(f)/2 min(py) 1.1*max(real(py(plotrange)))]), 
else
    axis([0 max(f) min(py) 1.1*max(real(py(plotrange)))]), 
end