Problem recovering the filtered signal using IFFT

4 次查看(过去 30 天)
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

回答(1 个)

Star Strider
Star Strider 2018-11-10
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.
  2 个评论
Artem Smirnov
Artem Smirnov 2018-11-10
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
Star Strider 2018-11-10
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Digital Filter Design 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by