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

42 views (last 30 days)
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)|')

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)|'); ##### 2 CommentsShowHide 1 older comment
Meg Noah on 9 Jan 2020
You are quite welcome!

### Categories

Find more on Discrete Fourier and Cosine Transforms in Help Center and File Exchange

### Tags

No tags entered yet.

### Community Treasure Hunt

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

Start Hunting!