Is my Fourier Transform code correct?

Hello.

I am using FFT to denoise a signal. However, my recovered signal does not match exactly my original signal.

Could someone check my code and let me know if I am doing something wrong?

clc;
close all;
clear all;
dt=0.001;
t=0:dt:1;

%CREATE AND PLOT INPUT SIGNAL 
signal=sin(2*pi*50*t)+sin(2*pi*120*t); % Input signal
figure('Name','1) Original Signal','NumberTitle','off'); %Create window to plot.
plot(t,signal); % Plot input "signal" as function of "t"
title('Original Signal');
ylabel('Amplitude');
xlabel('Time');
% ylim([-3 3]); %Set limit of y axis.

%CREATE AND PLOT NOISE SIGNAL
noise=2.5*randn(size(t)); % Noise
% subplot(3,1,2)
figure('Name','2) Noise','NumberTitle','off');
plot(t,noise);
title('Noise')
ylabel('Amplitude');
xlabel('Time');

%CREATE AND PLOT NOISY SIGNAL
noisysignal=signal+noise; % Noise signal as sum of input signal plus noise.
figure('Name','3) Noisy Signal','NumberTitle','off');
plot(t,noisysignal);
title('Noisy Signal');
ylabel('Amplitude');
xlabel('Time');

%CALCULATE AND PLOT MAGNITUDE (FOURIER TRANSFORM) OF NOISY SIGNAL
T=fft(noisysignal);  %Fourier Transform of "noisy signal"
fs = 1000; %Sample frequency
f = (0:length(T)-1)*fs/length(T); % Create frequnecy bins.
figure('Name','4) Amplitude Spectrum','NumberTitle','off');
plot(f,abs(T)); % Plot magnitude.
title('Amplitude Spectrum');
ylabel('Amplitude');
xlabel('Frequency (Hz)');

%CREATE AND PLOT CLEAN AMPLITUDE SPECTRUM
indices=abs(T)>300; % Variable "indices" only takes "abs(T)" values greater than 300
T=indices.*T; % Multiplied only 2 values or peaks (indices) by variable "T" to get a new Fourier Transform with only two values.
figure('Name','5) Clean Amplitude Spectrum','NumberTitle','off');
plot(f,abs(T)); % New plot of amplitude, will show only 2 amplitude values.
title('Clean Amplitude Spectrum');
ylabel('Amplitude');
xlabel('Frequency (Hz)');

%CALCULATE AND PLOT INVERSE FOURIER TRANSFORM.
ffilt=ifft(T); %Inverse Four Transform of only 2 amplitude values.
figure('Name','6) Recovered Signal','NumberTitle','off');
plot(t,ffilt);
title('Recovered Signal');
ylabel('Amplitude');
xlabel('Frequency (Hz)');

Answer 1

[Sulaymon Eshkabilov]

What you are doing is OK. Very small corrections/adjustments can be introduced for figure(4) and (5) to plot half spectrum excluding repeating parts:

L = length(t);
T=fft(noisysignal);  %Fourier Transform of "noisy signal"
TF=abs(T)/L;
TF1 = TF(1:L/2+1);
TF1(2:end-1) = 2*TF1(2:end-1);
Fs = 1/dt; %Sample frequency
f = Fs*(0:(L/2))/L;
...

Answer 2

[Paul]

Why whould this approach be expected to exactly recover the original signal? The DFT of the original signal is not zero everwhere except at the frequencies where indices == true. Also, the noise is contributing at those frequencies as well.

[Angel Lozada]

Paul.

Thanks for your answe.

May I know the value for dt? What shlould be?

Regards.