Extracting certain frequencies from FFT results using a window function

I am trying to come up with a code to extract certain frequencies from the fft result. I have 1 Hz data (259200 seconds for three days) and would like to convert the time series information into the frequency domain. The window size I would like to use is 2048 and the window function is hamming. The analysis window is shifted without any overlapping. The data at the frequency around 0.03 Hz is picked up from the FFT result and plotted against the time series data (x-axis) of 3 days (259200 seconds).

I know how to apply the window function on the whole length of the data but just not when you you have something like 2048 which is shorter than the length of the signal. I thought of using a loop function? Or probably buffer? Is this the right direction?

Also, I have no idea how to pick up the data from the fft results based on the frequency. How do I know the frequnecy I want (0.03 Hz) is in which bin? Meaning, if the frequency I want is within a bin then I would like the whole bin.

Answer 1

[Star Strider]

Frequency dokmain filtering is possible, however not ideal.

If you want to do that, the fftfilt function is the best option. It requires designing a FIR filter first, however that is relativelystraightforward using designfilt.

[KC]

Hi Strider

Thanks for the suggestions. This is what I have so far.

I have spent some time learning and compelling the code below. I have attached a txt file and if you could please point me in the right direction as I am afarid I haven't applied the window function or the fft correctly with respect to my question? I am I missing something?

data = readtable("data.txt");
t = [1:length(data.Var1)];

bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',length(data.Var1), 'Window','hamming');

y = fftfilt(bpfilt,data.Var1);    % filters the data in vector x with the filter described by coefficient vector b.

plot(t,abs(y))

Could you also please explain why frequency doamin is not ideal?

Thanks in advance for taking your time out to aid me in my understanding.

[Star Strider]

The filter appears to be working correctly.

My analysis —

data = readtable('data.txt');

bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',length(data.Var1), 'Window','hamming');

figure
freqz(bpfilt.Coefficients, 1, 2^16)
sgtitle('Filter Bode Plot')

% set(subplot(2,1,1), 'XScale','log')
% set(subplot(2,1,2), 'XScale','log')

s = data{:,1};
t = linspace(0, size(data,1)-1, size(data,1));
[FTs1,Fv] = FFT1(s,t);
% Fv(end)

figure
plot(Fv/Fv(end), mag2db(abs(FTs1)))
grid
Ax = gca;
% Ax.XScale = 'log';
xlabel('Normalised Frequency (rad/sample)')
ylabel('Magnitude (dB)')
title('Fourier Transform Of Original Signal')
% xlim([0 5.5E-5]*2*pi)
ylim([-200 50])

y = fftfilt(bpfilt,data.Var1);    % filters the data in vector x with the filter described by coefficient vector b.

[FTs1,Fv] = FFT1(y,t);
% Fv(end)

figure
plot(Fv/Fv(end), mag2db(abs(FTs1)))
grid
Ax = gca;
% Ax.XScale = 'log';
xlabel('Normalised Frequency (rad/sample)')
ylabel('Magnitude (dB)')
title('Fourier Transform Of Filtered Signal')
ylim([-200 50])

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

Plotting both the filtered and unfiltered signals on the same scale, however in different axes, seems to me to demonstrate that the filter is working as designed.

.

[KC]

Thank you Strider for the above. I do however, have a couple of questions if you would be so kind to answer them.

Firstly, in the designfilt part of the code, there wasn't any option to control the window size. I wanted to apply 1024 (2^11). Can the 2^16 be changed to 2^11?

Secondly, am I correct to understand that y = fftfilt(bpfilt,data.Var1) gives me the filtered results (between 0.027Hz and 0.033Hz) from the fft of the original signal? and this can then be plotted against the time?

Lastly, could you please explain why Frequency domain filtering is not ideal? Any links to resources will also be appreciated.

Thank you again for helping me out.

[Star Strider]

My pleasure!

Firstly, in the freqz call, the ‘2^16’ argument is the lenght of the fft the function calculates, and longer fft lengths result in greater frequency resolution. It has nothing to do with the window size of the filter.

Secondly, yes, and it can be realised by plotting ‘y’ as a function of the ‘t’ vector that I created (because my ‘FFT1’ function needs it). (Ideally, the time vector should be divided by the sampling frequency of the signal, however no sampling frequency was provided, so I just took the sampling frequency as 1 sample per ttime unit. The filters are all in terms of radian frequencies, so that works)

Lastly, frequency domain filtering is the rule in hardware filtering, however sampled signals exist in the time-domain, and time-domain filtering was developed for that purpose. The reason frequency domain filtering is not ideal for time-domain signals is that the signal and fflter need to be transformed back into the frequency domain to do the filtering. More computations create more inaccuracies in the eventual result, and of course, take more time. (A time-domain filter is essentially a frequency domain filter archetype transformed by taking its inverse Fourier transform, and then it is convolved with the time-domain recorded signal to do the filtering. In practice it is much more complicated than that, however that is the essential approach.)

[KC]

Hi Strider, thanks for the explanation above. It got me thinking a bit about time-domain filtering and I have a few more quesions regarding that with respect to frequency domain filtering.

Firstly, after we have extracted the frequency range (0.027-0.033 Hz) from the original signal, will the shape of the signal look similar to that of the original graph? Or is this a coincidence?

Secondly, let's say I did a time domain filtering with the butter function, should I also see a similar shape when compared to the fft filtered signal?

Lastly, I plotted all three graphs togther to see any differances or similarities as a result of filtering. Can I ask why there are differences even though we extracted the same frequency range?

Thank you again for sparing your time to help me understand.

data = readtable('data.txt');


%% FFT filter 

bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',length(data.Var1), 'Window','hamming');

s = data{:,1};
t = linspace(0, size(data,1)-1, size(data,1));
[FTs1,Fv] = FFT1(s,t);

y = fftfilt(bpfilt,data.Var1);    % filters the data in vector x with the filter described by coefficient vector b.

% plot(t,abs(y))

subplot(311)
plot(t,data.Var1)
title('Original Signal')

subplot(312)
plot(t,y)
title('FFT obtained Signal between 0.027Hz - 0.033Hz')


%% Butterworth filter 

input = data(:,1);

[k,c] = butter(1,[0.027 0.033],'bandpass'); % 2n so order 1

a=table2array(input);
a(isnan(a))=0;

y1 = filtfilt(k,c,a);

subplot(313)
plot(t,y1)
title('Butterworth Filtered Signal')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hamming(L), NFFT)/sum(hamming(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

[Star Strider]

My pleasdure!

Firstly, the extracted frequencies are so small in comparison to the frequencies in the signal, that the filter does not remove much. They essentially look the same.

Secondly, yes. The shapes will be similar. I did that experiment in the next-to-last plot.

Lastly. because the filter characteristics are definitely not the same, as can be seen in the freqz plots of each of them, and the last plot comparing the spectra of the filtered signals.

Filtering that narrow a range does not appear to make much difference in the filtered singal when compared to the unfiltered signal.

Testing that —

data = readtable('data.txt');


%% FFT filter 

bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',length(data.Var1), 'Window','hamming');

figure
freqz(bpfilt.Coefficients, 1, 2^18)

s = data{:,1};
t = linspace(0, size(data,1)-1, size(data,1));
[FTs1,Fv] = FFT1(s,t);

y = fftfilt(bpfilt,data.Var1);    % filters the data in vector x with the filter described by coefficient vector b.

% plot(t,abs(y))

figure
subplot(311)
plot(t,data.Var1)
title('Original Signal')

subplot(312)
plot(t,y)
title('FFT obtained Signal between 0.027Hz - 0.033Hz')


%% Butterworth filter 

input = data(:,1);

[k,c] = butter(1,[0.027 0.033],'bandpass'); % 2n so order 1



a=table2array(input);
a(isnan(a))=0;

y1 = filtfilt(k,c,a);

subplot(313)
plot(t,y1)
title('Butterworth Filtered Signal')

figure
freqz(k, c, 2^18)

figure
plot(t, a, 'DisplayName','Unfiltered')
hold on
plot(t, y, 'DisplayName','FFT Filtered')
plot(t, a-y, 'DisplayName','Difference')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

figure
plot(t, a, 'DisplayName','Unfiltered')
hold on
plot(t, y1, 'DisplayName','Butterworth Filtered')
plot(t, a-y1, 'DisplayName','Difference')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

figure
plot(t, y, 'DisplayName','FFT Filtered')
hold on
plot(t, y1, 'DisplayName','Butterworth Filtered')
plot(t, y-y1, 'DisplayName','Difference')
hold off
grid
xlabel('Time')
ylabel('Amplitude')
legend('Location','best')

[FTy,Fv] = FFT1(y,t);
[FTy1,Fv] = FFT1(y1,t);


figure
semilogx(Fv, mag2db(abs(FTy)*2), 'DisplayName','FFT Filtered')
hold on
plot(Fv, mag2db(abs(FTy1)*2), 'DisplayName','Butterworth Filtered')
plot(Fv, mag2db(abs(FTy))-mag2db(abs(FTy1)), 'DisplayName','Difference')
hold off
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Spectrum Comparison')
legend('Location','best')

% xlim([0 0.1])


function [FTs1,Fv] = FFT1(s,t)
s = s(:);
t = t(:);
L = numel(t);
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)).*hamming(L), NFFT)/sum(hamming(L));
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
FTs1 = FTs(Iv);
end

.

[Paul]

Can someone explain the SampleRate parameter in the call to designfilt?

data = readtable('data.txt');
%% FFT filter
bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',length(data.Var1), 'Window','hamming');

The frequency response of this filter is, assuming it's actually applied to a signal with sampling period Ts = 1 sec/sample (Fs = 1 Hz), is

[h1,f] = freqz(bpfilt.Coefficients,1,[],1);

Based on the question, it would seem that that the SampleRate parameter should be 1 sample/sec, or 1 Hz.

bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',1, 'Window','hamming');

With frequency response

h2 = freqz(bpfilt.Coefficients,1,f,1);

figure
subplot(211)
plot(f,db(abs(h1)),f,db(abs(h2))),grid
subplot(212)
plot(f,180/pi*angle(h1),f,180/pi*angle(h2)),grid
xlabel('f (Hz)')

[KC]

@Paul I think you maybe right since the documentation says "Sample Rate is the frequency at which the filter operates.". I think I may have overlooked it. Thanks for pointing it out.

[KC]

@Star Strider Thanks for all your help so far. I did had a couple of more questions and concerns if you could further aid me in my understanding.

Firstly, I think I overlooked on the “SampleRate” from above. We have to change it to 1 because the data is 1 Hz, right?

Secondly, I am trying to wrap my head around why the graph titled ‘FFT obtained Signal between 0.027Hz - 0.033Hz’ does not have similar amplitudes as that of the graph titled ‘Butterworth Filtered Signal’? How can the amplitudes be so far apart when we are looking at the same range of frequencies ? I get that these are different methods (frequency domain and time domain) but shouldn’t the amplitudes from the same range of frequencies be somewhat close together?

Thirdly, am I correct to understand that the ‘fftfilit’ used fft to transfer the original signal to the frequency domain and then used ifft to extract the frequency range between 0.027 - 0.033 Hz? The “y” we got is the amplitude of the fft?

Lastly, when applying this analysis to other lengths of data, such as two weeks or even a month, what I have stumbled upon is that ‘fftfilit’ gives nan values for them. Only when changing the ‘filter order’ in the ‘designfilit’ to 5 was I able to get rid of those nan’s. However, I wasn’t able to do that with the three month long data. Is there any workarounds to this? The ‘butter’ function, on the other hand, did not give any nan values. Do you think this is a better way?

Thanks again for helping me out.

[Star Strider]

My pleasure!

Firstly, I set the sampling frequency to 1 because the actual sampling frequency was not provided, and I like to have a time vector for the time domain plots. With careful coding, everything else can immediately adapt to whatever the sampling frequency actually is.

Secondly, the filter designs are different, and they have different characteristics, as can be seen in the freqz plots. That is not always as significant as it is here, with the pasband being as narrow as it is. Even if the freqz plot is zoomed significantly, it is difficult to see the passband in the freqz plot of the FIR filter —

data = readtable('data.txt');

%% FFT filter 

bpfilt = designfilt('bandpassfir', ...
    'FilterOrder',20,'CutoffFrequency1',0.027, ...
    'CutoffFrequency2',0.033,'SampleRate',length(data.Var1), 'Window','hamming');

figure
freqz(bpfilt.Coefficients, 1, 2^18, 1)
set(subplot(2,1,1), 'XLim',[0 0.25])
xline(subplot(2,1,1),[0.027 0.033]/2, '-r')
set(subplot(2,1,2), 'XLim',[0 0.25])
xline(subplot(2,1,2),[0.027 0.033]/2, '-r')

Although it is readily apparent in the Butterworth Bode plot. Staying in the frequency domain, even increasing the filter order significnatly does not make much difference in the FIR filter Bode plot.

Thirdly, essentially yes, although it is much more complicated than that. See the documentation section on Algorithms to understand how the fftfilt function works.

Lastly, I certainly have no objection to using FIR filters, and I routinely use them when I want to use multiple passbands or stopbands (‘comb’ filters) since that is more efficient than using IIR filters in series or parallel to do the same task. However for applications such as this, I routinely use IIR filters (specifically elliptic filters) for their computational efficiency. I would definitely not use a FIR filter for this application. I am also not certain what NaN values you are referring to. The NaN values are the result of 0/0, Inf/Inf, or NaN values in the data or an unstable filter (characteristically seen when using a transfer function realisation of an IIR filter, instead of a second-order-section realisation). Again, I would use a FIR filter for these data, and I recommend an elliptic design. The bandpass function can do that for you, simply by calling it as:

Fs = 1;

[y_filt, df] = bandpass(data.Var1, [0.027 0.033], Fs, 'ImpulseResponse','iir');

The ‘df’ output is a digital filter object that encodes an efficient elliptic filter. You can use it subsequently as:

y2_filt - filtfilt(df, y2);

to filter different signals with the same sampling frequency.

Again, my pleasure.

.

[KC]

@Star Strider thanks for the above. I just have a couple of more questions.

Firstly, is there a typo in the last part of your code above? More specifically, ' y2_filt - filtfilt(df, y2)'? If there isn't any mistake please could you elaborate more on the code.

Secondly, is there anyway we can control the window size in the fft frequency domain? I read there are a couple of ways to do that by multiplying the original signal to the hamming function? Is there a way to do this? I think maybe controlling the window size just may remove the NaN values after the doing the fft.

Thank you again for all your help. Much apprericated.

[Star Strider]

My pleasure.

The negative should be a n equal sign:

y2_filt = filtfilt(df, y2);

Windowing a Fourier transform corrects for the numeric Fourier transform being finite in length, as opposed to its actually being calculated over . The chosen window size (length) is the size (length) of the original time-domain vector.

The window function has no effect on the NaN values if they are in the data or if they are the result of filter instablility. It will not correct for them or eliminiate them.