how to smooth or filter the signal like this?
7 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
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!