How to use bandpass filter.

I have a data set which includes 2 columns: one is time (in month) and one is values that need analysics. I need to used bandpass filer to filter signal and obtain quasi-beninial oscillation (QBO) as figure 6(in attach). The period of this QBO is from 22 to 34 months, centered at 28 months. I try wirte a code but I don't belive on it. I would like everyone to check help me. My data, code, figure in attach. I hope everyone can help me.

Answer 1

[Mathieu NOE]

Hello

no bugs or special issues , but minor improvements

I prefer to use the IIR filters (like butter) used with filtfilt

clc;clear all;close all
load RTEC_N1.txt
time=RTEC_N1(:,1);
RTEC=RTEC_N1(:,2);

%% dt 
dt = mean(diff(time));
% fnorm_12=[1./28 1./24];
Fs=1/dt;
fc1 = 1./30;
fc2=1./22;

Wn = (2/Fs)*[fc1 fc2];
% Wn = [(2/Fs)*fc1 (2/Fs)*fc2];
% b = fir1(10,Wn,'low',kaiser(11,3));
% wn=[1./30 1./22];
% [b,a] = fir1(75,Wn);
% [b,a] = fir1(80,Wn,'bandpass');
% wn=1./(fs/2)/26
[b,a]=butter(2,Wn);
yy = filtfilt(b,a,RTEC);
% yy = filter(b,a,RTEC);
% figure(1);
% plot(time,RTEC,time,y,time,yy)
% legend('Original','filtfilt','filter');grid on;
plot(time,RTEC,'b',time,yy,'r');grid on;hold on;
fid1=fopen('RTIME_Sfiltertry1.txt','w');
for i=1:length(time)
    fprintf(fid1,'%10.3f\t%10.4f\r\n',time(i),yy(i));
end
fclose(fid1);

[Dung Nguyen]

Dear Mathieu NOE!

Thanks for your answer. I have run your code. It has not bug but the result seem to not reflect that I expect as figure 6 (in attach). In figure 6 the author used Fourier filtered. The period of this QBO is from 22 to 34 months, centered at 28 months (I don't know by which method put these periods into the filter). In my understanding, T1=24 months, T2=34 months corresponding fc1=1./34, fc2=1./22. Is it correct? We need using a filter which can keep periods or frequencies as above. I have attached the paper which used Fourier filtered to obtain figure 6. Could you show to help me? By another code using Fourier filter. I have struggled for a long time without finding the most optimal method. Thanks so much.

[Mathieu NOE]

hello

so when I am in doubt I do a fft analysis first to see where the signals has strong peaks in the spectrum .. here we have multiple peaks , and maybe the one around period = 26 months is the one you are looking for; then you can define the best filter parameters ...

updated code below :

clc;clear all;close all
load RTEC_N1.txt
time=RTEC_N1(:,1);  % time (in month)
RTEC=RTEC_N1(:,2);
%% dt 
dt = mean(diff(time));
Fs=1/dt;  % unit : 1 / month

%% first FFT analysis

    NFFT = length(RTEC);    % so freq resolution is highest
    OVERLAP = 0;

    [freq, sensor_spectrum] = myfft_peak(RTEC,Fs,NFFT,OVERLAP);
    % remove data for freq < 1/40 month and freq > 1/10 month ;
    ind=find(freq>1/40 & freq<1/10);
    freq = freq(ind);
    sensor_spectrum = sensor_spectrum(ind);
    %%%%
    period = 1./freq ; % period in month
    figure(1),plot(period,sensor_spectrum,'b');
    title('FFT Spectrum');
    xlabel('Frequency ( 1 / month)');ylabel('Amplitude');
    xlabel('Period ( month)');ylabel('Amplitude');


%%% then filtering
fc = 1/26;  % estimated central frequency 
fc1 = fc/1.2; % lower  and higher cut off filter frequencies are symetrical from fc with a factor of 1.2
fc2 = fc*1.2;
Wn = (2/Fs)*[fc1 fc2];

[b,a]=butter(2,Wn);
yy = filtfilt(b,a,RTEC);

figure(2),plot(time,RTEC,'b',time,yy,'r');grid on;hold on;
% write results
writematrix([time(:) yy(:)], 'RTIME_Sfiltertry1.txt',"Delimiter","tab");


%% second FFT analysis

    NFFT = length(RTEC);    % so freq resolution is highest
    OVERLAP = 0;

    [freq, sensor_spectrum1] = myfft_peak(RTEC,Fs,NFFT,OVERLAP);
    [freq, sensor_spectrum2] = myfft_peak(yy,Fs,NFFT,OVERLAP);

    %%%%
    period = 1./freq ; % period in month
    figure(3),semilogx(period,sensor_spectrum1,'b',period,sensor_spectrum2,'r');
    title('FFT Spectrum');
    xlabel('Frequency ( 1 / month)');ylabel('Amplitude');
    legen('raw data','filtered data');
    xlabel('Period ( month)');ylabel('Amplitude');
    
    

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[Dung Nguyen]

Hello

Thanks for your code. It is useful for me. I still have a problem and need your help. I want to use wavelet transform to find quasi- biennial period and I have used a code which I downloaded, but the result is not good. Because I can't see a quasi- biennial period in the result. Would you check to help me? Thanks so much.

[Mathieu NOE]

hello

I have to admit that both the regular fft and the wavelet analysis do not really highlight any important contribution of a 28 period signal in the raw data... I doubt that the fft or the wavelet tool are in cause , so what does you feel so sure that the raw data does contain that 28 month period cycle ??

FFT analysis show more dominant peaks than a 28 month period

Wavelet idem

raw data : so irregular , one would be cautious about a quasi- biennial period signal (??)

[Mathieu NOE]

last figure missing :

[Dung Nguyen]

Hello

Thanks for your replying. I have tried to find period by Lomb- Scarge periodogram (in attached). I think the same you that 26 or 28 month period is weaker than other period. "so what does you feel so sure that the raw data does contain that 28 month period cycle ??"....> I watched the temporal variation by my eyes in raw data. For my eyes, it looks to have a weak quasi-bianual variation. May be it is weak so we difficultly find out it in analysis as wavelet. It is true that I can't do what more. Do you have other ideal? Hope that You can show to help me. Thanks so much.

[Mathieu NOE]

hello

well, your eyes are more powerful than any spectral analysis tool !! some tools are more appropriate than others for non stationnary signals (like wavelet vs fft) , but as already said above, there is some trace of the 26-28 month period signal but it's weak; does it means there is a problem with the data ? honestly , i cannot say because my field of expertise is more industrial vibration analysis and control , rather than weather or oceanographic data .

I may come to the point where my help is coming to and end - sorry for that

[Dung Nguyen]

Hello

Thanks for your help.