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

43 次查看(过去 30 天)
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
Meg Noah 2020-1-9
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 个评论
Tamura Kentai
Tamura Kentai 2020-1-9
Yes, from the output waveform, I think it's correct.
Thank you !!!

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Digital Filter Analysis 的更多信息

标签

尚未输入任何标签。

产品

Community Treasure Hunt

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

Start Hunting!

Translated by