alternativeredrawfourierfilter.m

来自「A fourier filter for time-series signals」· M 代码 · 共 53 行

function ry=RedrawFourierFilter(xvector,yvector,centerfrequency,filterwidth,filtershape,mode)
% Alternative version of RedrawFourierFilter, with original signal (blue)
% and filtered signal (red) displayed superimposed in the top window.
% Computes and plots fourier filter for signal yvector.  centerfrequency
% and filterwidth are the center frequency and width of the
% pass band. 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.
% To display original signal and filtered signal in the same window, 
% rename AlternativeRedrawFourierFilter.m to RedrawFourierFilter.m 
%  T. C. O'Haver ([e-mail]),  version 1.4, October 29, 2006
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 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 Transfork of signal
ry=real(ifft(ffy));
% Plot original and filtered signal in top plot, power spectrum and filter
% in lower plot
subplot(2,1,1)
plot(xvector,yvector,xvector,ry,'r')
title(['BLUE = ' num2str(signalstring) '        RED = Filtered signal'])

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