Problem recovering the filtered signal using IFFT

Hello! I have geomagnetic field values, 1 value each month, for several years. What I want to do with the data is to erase yearly-signal and smaller period signals (frequency 0.0027 Hz and everything above that)from the initial data using the Butterworth filter in the frequency domain, and then restore the filtered signal using the IFFT. However, I get very strange output when I do it.

The spectrum seems to be correct, but how can I recover the signal after filtering it? I attach the file 'H.mat', my code and the screenshot.

What can be the problem?

load('H.mat');
time=H(:,1);
field=H(:,2);

n = length(time);
Ts = 1/30*24*60;                  % Sampling Interval (1 month)
Fs = 1/Ts;                        % Sampling Frequency (1 value per month)

Fn = Fs/2;                         % Nyquist Frequency
nfft = 2^nextpow2(n);
ft = fft(field,nfft)/n;
Fv = linspace(0, 1, fix(n/2)+1)*Fn;     % Frequency Vector
Iv = 1:length(Fv);                      % Index Vector

figure(1)
subplot(3,1,1)
plot(time, field,'k.-')
xticks(2000:2016)
title('Initial field')
grid

subplot(3,1,2)
plot(Fv, abs(ft(Iv))*2)
title('Fourier spectrum. I need to filter out f=0.002 and higher ones')
grid

%Now the filter
[b,a]=butter(6,0.0027);
filtered_field=filter(b,a,ft);
%recover the field in time-domain:
field_output=ifft(filtered_field);

subplot(3,1,3)
plot(field_output)
title('Field after IFFT ?... ')
grid on

Answer 1

[Star Strider]

Your filtering approach is incorrect.

Try this:

load('H.mat');
time=H(:,1);
field=H(:,2);

n = length(time);

Ts = mean(diff(time));
Fs = 1/Ts;
Fn = Fs/2;

nfft = 2^nextpow2(n);
ft = fft(field,nfft)/n;
Fv = linspace(0, 1, fix(n/2)+1)*Fn;     % Frequency Vector
Iv = 1:length(Fv);                      % Index Vector

figure(1)
subplot(3,1,1)
plot(time, field,'k.-')
xticks(2000:2016)
title('Initial field')
grid

subplot(3,1,2)
plot(Fv, abs(ft(Iv))*2)
title('Fourier spectrum. I need to filter out f=0.002 and higher ones')
grid

Wp = 0.0020/Fn;                                             % Passband Frequency Vector (Normalised)
Ws = 0.0019/Fn;                                             % Stopband Frequency Vector (Normalised)
Rp =   1;                                                   % Passband Ripple (dB)
Rs =  50;                                                   % Stopband Attenuation (dB)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                             % Calculate Filter Order
[z,p,k] = ellip(n,Rp,Rs,Wp);                                % Default Here Is A Lowpass Filter
[sos,g] = zp2sos(z,p,k);                                    % Use Second-Order-Section Implementation For Stability
field_output = filtfilt(sos,g,field);                       % Filter Signal (Here: ‘field’)

subplot(3,1,3)
plot(time, field_output)
title('Field Filtered')
grid on

figure(2)
freqz(sos, 2^14, Fs)                                        % Bode Plot Of Filter
% set(subplot(2,1,1), 'XLim',[0 15])                          % Optional, Change Limits As Necessary
% set(subplot(2,1,2), 'XLim',[0 15])                          % Optional, Change Limits As Necessary

You may need to experiment with the filter passband if you want a different result. That should not be difficult.

[Artem Smirnov]

Thanks! But the result is the exact opposite of what I should get... When I applied it, I got the perfect yearly signal, while the main aim of the filtering is to delete exactly this signal... Do you know if there is any lowpass filter that can do it?

[Star Strider]

My pleasure.

A lowpass filter passes only the frequencies below the cutoff frequency, here 0.002, so the output will be from 0 to 0.002. A highpass filter will pass everything above the cutoff frequency, here from 0.002 to the Nyquist frequency.

To change the filter to a highpass design, change only the ellip call to:

[z,p,k] = ellip(n,Rp,Rs,Wp,'high');                         % Default Here Is A Lowpass Filter

That worked when I tried it.