About FFT of sin(x) / x , how to plot ?

79 views (last 30 days)
Tamura Kentai
Tamura Kentai on 9 Jan 2020
Commented: Meg Noah on 9 Jan 2020
Hi :
I can't find solutions after searched.
I tried to plot a FFT spectrum of the sinc function, y(x) = sin(x) /x, but I can't get the correct output, the 'FFT'ed value as Y here, I saw the content of 'Y' is almost NaN.
Can anyone help solve this ?
I have tried to comment the line 'S = Sn./t;' ,and replaced it with the next line 'S = Sn;', it successfully ran out the FFT spectrum with a obvious component at F (x-axis) = 100 (Hz). I guess the most code in it works fine except the math expression "S = sin(t) / t" here , it's " Sn = sin(2*pi*1e2*t); S = Sn ./t; ".
I tried to plot the frequency response of 'brick wall'.
Thank you very much.
% show the " y=sin(t) / t " in FFT spectrum
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 15e4; % Length of signal
% t = (0:L-1)*T; % Time vector
t = (-L/2:L/2-1)*T; % Time vector
% Form a signal containing ' y=sin(t) / t ' ,here 'S' is the Signal as 'y' in the math expression.
Sn = sin(2*pi*1e2*t);
S = Sn ./t;
% S = Sn;
% Corrupt the signal with zero-mean white noise with a variance of 4.
% X = S + 2*randn(size(t));
X = S;
% Plot the noisy signal in the time domain. It is difficult to identify the frequency components by looking at the signal X(t).
plot(1000*t(1:50),X(1:50))
title('Signal X(t) = sin(t) / t')
xlabel('t (milliseconds)')
ylabel('X(t)')
% Compute the Fourier transform of the signal.
Y = fft(X);
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% Define the frequency domain f and plot the single-sided amplitude spectrum P1. The amplitudes are not exactly at 0.7 and 1, as expected, because of the added noise. On average, longer signals produce better frequency approximations.
f = Fs*(0:(L/2))/L;
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Accepted Answer

Meg Noah
Meg Noah on 9 Jan 2020
Hi, I was able to get a plot by removing the discontinuity in the signal. Is this what you were expecting?
% show the " y=sin(t) / t " in FFT spectrum
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 15e4; % Length of signal
% t = (0:L-1)*T; % Time vector
t = (-L/2:L/2-1)*T; % Time vector
% Form a signal containing ' y=sin(t) / t ' ,here 'S' is the Signal as 'y' in the math expression.
Sn = sin(2*pi*1e2*t);
S = Sn ./t;
% remove discontinuity
idx = find(isnan(S));
S(idx) = 1;
% S = Sn;
% Corrupt the signal with zero-mean white noise with a variance of 4.
% X = S + 2*randn(size(t));
X = S;
% Plot the noisy signal in the time domain. It is difficult to identify the frequency components by looking at the signal X(t).
figure()
subplot(3,1,1);
plot(t,X);
xlim([-0.2 0.2]);
title('Signal X(t) = sin(t) / t with zero-mean white noise \sigma=4');
xlabel('t (s)')
ylabel('X(t)');
subplot(3,1,2);
plot(1000*t(1:50),X(1:50))
title('Signal X(t) = sin(t) / t with zero-mean white noise \sigma=4');
xlabel('t (milliseconds)')
ylabel('X(t)');
% Compute the Fourier transform of the signal.
Y = fft(X);
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% Define the frequency domain f and plot the single-sided amplitude spectrum P1. The amplitudes are not exactly at 0.7 and 1, as expected, because of the added noise. On average, longer signals produce better frequency approximations.
f = Fs*(0:(L/2))/L;
subplot(3,1,3)
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|');
Discontinuity.png
Explanation in this Youtube video:
  2 Comments
Meg Noah
Meg Noah on 9 Jan 2020
You are quite welcome!

Sign in to comment.

More Answers (0)

Tags

No tags entered yet.

Products

Community Treasure Hunt

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

Start Hunting!

Translated by