how to smooth or filter the signal like this?
3 views (last 30 days)
Show older comments
How to smooth the signal like this?
I want to do fourier transform of the data attached. I only need the sharp peak in the data, I don't know how to
filter the noise in the data. what about Hanning window or something else?
Your help would be highly appreciated.
Best regards
load('data.mat')
plot(t,dip)
0 Comments
Answers (4)
Bora Eryilmaz
on 19 Dec 2022
Edited: Bora Eryilmaz
on 19 Dec 2022
You can apply a low-pass filter to remove the noise in the signal before taking the Fourier Transform of the data:
load('data.mat')
plot(t,dip)
hold on
Fs = 1/(t(2)-t(1));
[b,a] = butter(4, 0.03/(Fs/2), 'low');
y = filter(b,a, real(dip));
plot(t,y)
By the way, your original data, dip, seems to have imaginary values in it. You would probably want to get rid of those.
0 Comments
Star Strider
on 19 Dec 2022
Edited: Star Strider
on 19 Dec 2022
I saved that code.
I did not save the non-code text that went with it. I am not certain that I can reproduce the image that you posted, because I do not understand what it is or how it was calculated. I plotted the Fourier transform three times here, once with the absolute amplitude, once with the amplitude in decibels, and once with the log amplitude.
Reposted —
LD = load(websave('data','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1236142/data.mat'));
t = LD.t(:);
dip = LD.dip(:);
figure
plot(t, real(dip))
grid
xlabel('t')
ylabel('AMplitude')
title('Original Signal')
xlim([min(t) max(t)])
Ts = mean(diff(t)); % Sampling Interval
Fs = 1/Ts; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
L = numel(t);
NFFT = 2^nextpow2(L); % For Efficiency
FTdip = fft(dip-mean(dip).*hamming(L),NFFT)/L; % Windowed Fourier Transform
Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
figure
plot(Fv, abs(FTdip(Iv))*2, 'DisplayName','Data')
grid
xlabel('Frequency')
ylabel('Magnitude (Absolute)')
title('One-Sided Fourier Transform')
% xlim([0 0.5])
[maxPeak,idx] = max(abs(FTdip(Iv))*2)
maxPeakFreq = Fv(idx)
passband = maxPeakFreq*[0.75 1.25]
xline(passband, '-r', 'DisplayName','Passband')
legend('Location','best')
figure
plot(Fv, mag2db(abs(FTdip(Iv))*2), 'DisplayName','Data')
grid
xlabel('Frequency')
ylabel('Magnitude (dB)')
title('One-Sided Fourier Transform')
% xlim([0 0.5])
xline(passband, '-r', 'DisplayName','Passband')
legend('Location','best')
figure
semilogy(Fv, abs(FTdip(Iv))*2, 'DisplayName','Data')
grid
xlabel('Frequency')
ylabel('Magnitude')
title('One-Sided Fourier Transform')
% xlim([0 0.5])
xline(passband, '-r', 'DisplayName','Passband')
legend('Location','best')
filtered_dip = bandpass(real(dip), passband, Fs, 'ImpulseResponse','iir');
figure
plot(t, filtered_dip)
grid
xlabel('t')
ylabel('Amplitude')
title('Filtered Signal')
xlim([min(t) max(t)])
This code selects the maximum frequency and defines a passband based on it. It then filters the signal using that passband. (The passband is also depicted in the frequency domain plot to illustrate what that filter passes.) I do not suggest a narrower passband, however shifting the existing passband or enlarging it (or both) are straightforward.
EDIT — (19 Dec 2022 at 20:36)
Looking at thid further, I am not certain how to get the result you want from this isgnal. The Savitzky-Golay filter (sgolayfilt) may be the way to go on it (another option being a FIR comb filter), however it is not obvious to me how to get the result depicted in your latest plot image.
% new_dip = lowpass(real(dip), 0.1, Fs, 'ImpulseResponse','iir')
sgfilt_dip = sgolayfilt(real(dip), 3, 127);
figure
plot(t, sgfilt_dip)
grid
xlabel('t')
ylabel('Amplitude')
title('Savitzky-Golay Filtered Signal')
xlim([min(t) max(t)])
FTdip = fft(sgfilt_dip-mean(sgfilt_dip).*hamming(L),NFFT)/L; % Windowed Fourier Transform
figure
semilogy(Fv, abs(FTdip(Iv))*2, 'DisplayName','Data')
grid
xlabel('Frequency')
ylabel('Magnitude')
title('One-Sided Fourier Transform Of Savitzky-Golay Filtered Signal')
% xlim([0 0.5])
xline(passband, '-r', 'DisplayName','Passband')
legend('Location','best')
Experiment with the sgolayfilt function. Also, it would help if you described in some detail the result you want, since that is not obvious to me.
.
0 Comments
Image Analyst
on 19 Dec 2022
Not sure what you consider a sharp peak. Is it anything above 0.25 or below -0.25? If so, you can just zero out everything with a magnitude less than 0.25:
load('data.mat')
dip(abs(dip) < 0.25) = 0;
plot(t, dip)
Then do whatever you want, such as filtering in the spectral domain.
0 Comments
See Also
Categories
Find more on Digital Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!