Is my Fourier Transform code correct?

1 view (last 30 days)
Angel Lozada
Angel Lozada on 15 May 2021
Commented: Angel Lozada on 18 May 2021
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)');

Answers (2)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 15 May 2021
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;
...

Paul
Paul on 15 May 2021
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.
  1 Comment
Angel Lozada
Angel Lozada on 18 May 2021
Paul.
Thanks for your answe.
May I know the value for dt? What shlould be?
Regards.

Sign in to comment.

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!