MATLAB Answers

Denoising the sinusoidal signal

7 views (last 30 days)
adjocket
adjocket on 20 May 2021
Commented: adjocket on 20 May 2021
How can I put the xn signal into the dee filter and draw it?
numSamples = 250;
t = linspace(0, 50, numSamples);
x = sin(2*pi*100*t) + 0.7*cos(2*pi*500*t);
xn = x + 5.*randn(size(t));
originalSpectrum = fft(x);
noisySpectrum = fft(xn);
subplot(4, 1, 1);
plot(abs(originalSpectrum), 'b-', 'LineWidth', 2);
grid on;
xlabel('t');
ylabel('x');
subplot(4, 1, 2);
plot(abs(noisySpectrum), 'b-', 'LineWidth', 2);
dee = designfilt('lowpassfir', 'PassbandFrequency', 500, ...
'StopbandFrequency', 600, ...
'PassbandRipple', 1, 'StopbandAttenuation', 100, ...
'SampleRate', 5000);

Answers (2)

Jonas
Jonas on 20 May 2021
use
dataOut = filter(dee,xn)
  2 Comments
Jonas
Jonas on 20 May 2021
sure? if i use your code and add my line it just works
edit: i tried that in matlab online. which matlab do you use? maybe the lowpass() function is also suitable for you

Sign in to comment.


Mathieu NOE
Mathieu NOE on 20 May 2021
hello
see my suggestion below
notice that your code was wrong regarding sampling frequency (barely at 5 Hz) so this is not coherent if you want to generate signals at 100 and 500 Hz (Shannon's theorem implies you must at least sample at twice the highest frequency => Fs > 1000 Hz)
also , usually in fft analysis, the spectrum amplitude are displayed in log scale or , which is equivalent , in dB scale - as here ; then it's much easier to see low and high values together , if you keep a linear y axis, the very small amplitudes remain invisible
the 0 dB level correspond to a sinewave of amplitude = 1
a sine wave with amplitude 0.7 would correspond to -3 dB in dB scale
a sine wave with amplitude 0.1 would correspond to -20 dB in dB scale
etc...
hope this helps
numSamples = 10000;
Fs = 2000;
dt = 1/Fs;
t = (0:numSamples-1)*dt;
x = sin(2*pi*100*t) + 0.7*cos(2*pi*500*t);
xn = x + 0.05.*randn(size(t));
% FFT Analysis
nfft = 1000;
Overlap = 0.5;
[freq_vector,originalSpectrum] = myfft_peak(x(:), Fs, nfft, Overlap);
[freq_vector,noisySpectrum] = myfft_peak(xn(:), Fs, nfft, Overlap);
dee = designfilt('lowpassfir', 'PassbandFrequency', 500, ...
'StopbandFrequency', 600, ...
'PassbandRipple', 1, 'StopbandAttenuation', 100, ...
'SampleRate', 5000);
xnf = filter(dee,xn);
[freq_vector,FilteredSpectrum] = myfft_peak(xnf(:), Fs, nfft, Overlap);
plot(freq_vector,20*log10(originalSpectrum), 'b-',...
freq_vector,20*log10(noisySpectrum), 'r--',...
freq_vector,20*log10(FilteredSpectrum), 'g-.', 'LineWidth', 2);
grid on;
xlabel('Frequency (Hz)');
ylabel('Amplitude (dB scale)');
legend('original Spectrum','noisy Spectrum','Filtered Spectrum');
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
s_tmp((1:samples),:) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
cor_coef = length(window)/sum(window);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_spectrum = 0;
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
fft_spectrum = fft_spectrum + (abs(fft(sw))*2*cor_coef/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only
end
fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end
  1 Comment
adjocket
adjocket on 20 May 2021
Actually this is the question. I'm trying to do as much as I can, but I couldn't figure it out. Is there a simpler way of making it?

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!