Plotting DTFT in Matlab

172 次查看（过去 30 天）
Tomasz Nowopolski
Tomasz Nowopolski 2022-11-26
编辑： Paul 2024-10-27，13:38
I'm trying to create a DTFT function in Matlab for my assignment.
I'm using the built-in FFT function to compare my results and wager if my results are correct (I'm assuming the results should be the same just with different error margins), but the result I'm getting is entirely different and I'm not sure why.
Main code:
L = 1000;
fg = 1000;
fs = 100000;
T = 1/fs;
tmin = 0;
K = 1;
n = tmin:L;
t = n*T;
%Base function calculation
x = xdp1(K, fg, t);
%Fourier transform calculation
w = t;
%X = dtft(x, t);
%X = dtft2(x, w, n);
X = dtft3(x, w, n);
X_FFT = fft(x);
% Top plot
tiledlayout(3,1)
nexttile
plot(t, x)
ylim([min(x) max(x)])
xlim([0 max(t)])
title('Signal')
% Middle plot
nexttile
plot(t, abs(X));
title('DTFT')
% Bottom plot
nexttile
plot(t, abs(X_FFT))
ylim([min(abs(X_FFT)) max(abs(X_FFT))])
xlim([0 max(t)])
title('FFT')
Signal calc code (xdp1):
function x = xdp1(K, fg, t)
x = K * 2 * fg * sinc(2 * pi * fg * t);
end
DTFT calc code (dtft3):
function X = dtft3(x, w, n)
X = exp(-1i*w'*n) * x.';
end
My other attempts that minimized the Matlab exclusive syntax:
DTFT calc code (dtft2):
function X = dtft2(x, w, n)
for i=1:length(w)
X(i)=sum(x.*exp(-1i*w(i)*n));
end
end
DTFT calc code (dtft):
function X = dtft(x, w)
X = zeros(1, length(x));
for n=1:length(x)
X = X + x(n) * exp(-1i * w * n);
end
end
Result plot:

请先登录，再回答此问题。

回答（2 个）

Ishan
Ishan 2022-11-29
Hi Tomasz,
I understand you are trying to cross verify your DTFT function. You can refer to the functions in the file exchange attached below
Digital Signal Processing -- Discrete-time Fourier Transform (DTFT) - File Exchange - MATLAB Central (mathworks.com)
Hope this helps out!
Paul
Paul 2024-10-26，20:53
编辑：Paul 2024-10-27，13:38
Ran in:
L = 1000;
fg = 1000;
fs = 100000;
T = 1/fs;
tmin = 0;
K = 1;
n = tmin:L;
t = n*T;
%Base function calculation
x = xdp1(K, fg, t);
Based on how the dtft* functions are defined, the discrete-time frequency vector should be
%Fourier transform calculation
% w = t;
w = 2*pi*n/numel(n); % rad/sample
The function dtft() should be called with w, not t. The function dtft() is corrected below. Test all three variations
%X1 = dtft(x, t);
X1 = dtft(x, w);
X2 = dtft2(x, w, n);
X3 = dtft3(x, w, n);
The DTFT for a finite duration signal can also be computed with freqz
X4 = freqz(x,1,w);
X_FFT = fft(x);
figure
% Top plot
tiledlayout(3,1)
nexttile
plot(t, x)
ylim([min(x) max(x)])
xlim([0 max(t)])
xlabel('t')
title('Signal')
% Middle plot
nexttile
%plot(t, abs(X));
plot(w,abs(X1))
hold on
plot(w,abs(X2))
plot(w,abs(X3))
plot(w,abs(X4))
xlabel('\omega (rad/sample)')
ylabel('Amplitude')
title('DTFT')
% Bottom plot
nexttile
Plot vs w, not t
%plot(t, abs(X_FFT))
plot(w, abs(X_FFT))
ylim([min(abs(X_FFT)) max(abs(X_FFT))])
%xlim([0 max(t)])
xlabel('\omega (rad/sample)');
ylabel('Amplitude')
title('FFT')
figure
% Top plot
tiledlayout(3,1)
nexttile
plot(t, x)
ylim([min(x) max(x)])
xlim([0 max(t)])
xlabel('t')
title('Signal')
% Middle plot
nexttile
%plot(t, abs(X));
plot(w,angle(X1))
hold on
plot(w,angle(X2))
plot(w,angle(X3))
plot(w,angle(X4))
xlabel('\omega (rad/sample)')
ylabel('Phase (rad)')
title('DTFT')
% Bottom plot
nexttile
%plot(t, abs(X_FFT))
plot(w, angle(X_FFT))
%ylim([min(abs(X_FFT)) max(abs(X_FFT))])
%xlim([0 max(t)])
xlabel('\omega (rad/sample)');
ylabel('Phase (rad)')
title('FFT')
%Signal calc code (xdp1):
function x = xdp1(K, fg, t)
x = K * 2 * fg * sinc(2 * pi * fg * t);
end
%DTFT calc code (dtft3):
function X = dtft3(x, w, n)
X = exp(-1i*w'*n) * x.';
end
%My other attempts that minimized the Matlab exclusive syntax:
%DTFT calc code (dtft2):
function X = dtft2(x, w, n)
for i=1:length(w)
X(i) = sum(x.*exp(-1i*w(i)*n));
end
end
%DTFT calc code (dtft):
function X = dtft(x, w)
X = zeros(1, length(x));
% x(1) corresponds to x[n = 0], so the sum over n starts at n = 0.
%for n=1:length(x)
for n = 0:numel(x)-1
X = X + x(n+1) * exp(-1i * w * n);
end
end

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by