Denoising the sinusoidal signal

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);

Answer 1

[Jonas]

use

dataOut = filter(dee,xn)

[adjocket]

gives this error (Function 'subsindex' is not defined for values of class 'digitalFilter')

[Jonas]

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

Answer 2

[Mathieu NOE]

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

[adjocket]

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?