DTFT possible on Matlab?
200 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hello everyone, I understand the usage of DFT but I would like to specifically perform a DTFT on a signal. Is it possible to do so in Matlab?
Thank you very much in advance!
0 Kommentare
Antworten (7)
Paul
am 5 Jun. 2023
Bearbeitet: Paul
am 7 Jun. 2023
Computing the DTFT of a signal in Matlab depends on
a) if the signal is finite duration or infinite duration
b) do we want the numerical computation of the DTFT or a closed form expression.
In the examples that follow, u[n] is the discrete time unit step function, i.e.,
u[n] = 1, n >= 0
u[n] = 0, n < 0
Case 1: finite duration, numerical computation
Let:
x1[n] = (0.5^n)*(u[n] - u[n - 8])
x2[n] = (1.5^n)*(u[-n] - u[-n + 5])
x[n] = x1[n] + x2[n];
x[n] is finite duration, i.e., x[n] = 0 for n < -4 and n > 7.
u = @(n) double(n>=0); % discrete unit step
x1 = @(n) (0.5.^n).*(u(n) - u(n-8));
x2 = @(n) (1.5.^n).*(u(-n) - u(-n - 5));
x = @(n) x1(n) + x2(n);
Numerical computation of the DTFT can be accomplished in a couple of ways. Perhaps the easiest is to use the Signal Processing Toolbox function freqz. However, using freqz in this way assumes that the first element of the first input to freqz corresponds to n = 0. But, in this example the first nonzero element of x[n] is at n = -4. We have to use the shift property of the DTFT after calling freqz to ensure we get the correct phase.
Let y[n] be a version of x[n] shifted to the right by m samples such that all of the non-zero elements of y[n] are at times n >= 0.
y[n] = x[n - m], m > 0
In our case, we have to shift x[n] to the right by at least 4 samples.
The DTFT of y[n] is related to the DTFT of x[n] as
Y(w) = X(w)*exp(-1j*w*m) (w in rad)
Therefore, we can use freqz to compute Y(w), then invert the above to compute X(w)
X(w) = Y(w)*exp(1j*w*m)
Compute x[n] over its duration
n = -4:7;
xval = x(n);
Compute freqz of xval, which are really the values of yval as far as freqz is concerned. Use the 'whole' option to evaluate around the entire unit circle.
[h,w] = freqz(xval,1,2048,'whole');
Phase adjustment for the time domain shift implicit to freqz
h = h.*exp(1j*w*4);
h is basically one period of the DTFT of x[n]. Plot it
figure
plot(subplot(211),w,abs(h)); grid, ylabel('Magnitude'),ylim([0 5]),xlim([0 2*pi])
plot(subplot(212),w,angle(h)); grid, ylabel('Phase (rad)'), xlabel('w (rad)'),xlim([0 2*pi])
Another option, which might be slightly more numerically robust, is to use fft with lots of zero padding to compute one period of the DTFT at lots of frequencies. fft assumes the first point in the sequence corresponds to n = 0, so we have to circshift after zero padding and then apply fft
h = fft(circshift([xval zeros(1,2048-numel(xval))],-4),2048);
w = (0:2047)/2048*2*pi;
figure
plot(subplot(211),w,abs(h)); grid, ylabel('Magnitude'),ylim([0 5]),xlim([0 2*pi])
plot(subplot(212),w,angle(h)); grid, ylabel('Phase (rad)'), xlabel('w (rad)'), xlim([0 2*pi])
Another option would be write a function to compute the DTFT sum directly from its definition at frequencies of interest.
Case 2: Finite duration, closed form expression.
This case is straightforward using the Symbolic Math Toolbox
clear all
syms n integer
u(n) = kroneckerDelta(n)/2 + heaviside(n); % discrete time unit step using default sympref
x1(n) = (sym(0.5)^n)*(u(n) - u(n-8));
x2(n) = (sym(1.5)^n)*(u(-n) - u(-n - 5));
x(n) = x1(n) + x2(n);
We know that x[n] is zero outside the interval -4 <= n <= 7, so we can use the DTFT sum over that interval explicitly
syms w
nval = -4:7;
X(w) = sum(x(nval).*exp(-1j*w*nval))
Plot one period
figure
fplot(subplot(211),abs(X(w)),[0 2*pi]), grid, ylabel('Magnitude'), ylim([0 5])
fplot(subplot(212),angle(X(w)),[0 2*pi]), grid, ylabel('Phase (rad)'), xlabel('w (rad)')
Case 3: Infinite duration, closed form expression
Let x[n] be similar to above, but instead of cutting x[n] off at n = -5 and n = 8, we'll let it run out to abs(n) = inf.
clear all
syms n integer
u(n) = kroneckerDelta(n)/2 + heaviside(n); % discrete time unit step using default sympref
x1(n) = (sym(0.5)^n)*(u(n));
x2(n) = (sym(1.5)^n)*(u(-n));
x(n) = x1(n) + x2(n);
Matlab's Symbolic Math Toolbox doesn't offer a symbolic DTFT function, but it does offer the unilateral z-transform via ztrans/iztrans. We can use ztrans to compute the bilateral z-transform of x1[n], which is the same as its unilateral transform because x1[n] is causal. We can use ztrans to to compute the bilateral z-transform of x2[n] using the z-transform time reversal property because x2[n] is strictly non-causal. Because the bilateral z-transforms of x1[n] and x2[n] both have a region of convergence that includes the unit circle, so does their sum, in which case the DTFT of x[n] can be computed by evaluating its bilateral z-tansform around the unit circle.
syms z
X1(z) = ztrans(x1(n),n,z); % z-transform of causal signal, ROC: abs(z) > 0.5
X2(z) = simplify(subs(ztrans(x2(-n),n,z),z,1/z)); % z-transform of non-causal signal, ROC: abs(z) < 1.5
X(z) = X1(z) + X2(z); % ROC: 0.5 < abs(z) < 1.5
syms w
X(w) = X(exp(1j*w)); % DTFT of x[n]
Plot one period
figure
fplot(subplot(211),abs(X(w)),[0 2*pi]), grid, ylabel('Magnitude'), ylim([0 5])
fplot(subplot(212),angle(X(w)),[0 2*pi]), grid, ylabel('Phase (rad)'), xlabel('w (rad)'), ylim([-0.2 0.2])
Case 4: Infinite duration, numerical computation
Because x[n] is infinite duration we can't really do a numerical computation, which would need an infinite number of samples. But x[n] decays to zero in both directions, so we can pick a duration of nlower <= n <= nupper, where for values of n outside those limits it's safe to assume that x[n] can be approximated as x[n] = 0.
We'll use upper and lower bounds of 20
nval = -20:20;
xval = double(x(nval));
Now we can use any of the numerical techniques from Case 1 with xval being elements of the finite duration approximation to x[n]
[hval,wval] = freqz(xval,1,2048,'whole');
hval = hval.*exp(1j*wval*20);
Plot one period
figure
plot(subplot(211),wval,abs(hval)); grid, ylabel('Magnitude'),ylim([0 5]),xlim([0 2*pi])
plot(subplot(212),wval,angle(hval)); grid, ylabel('Phase (rad)'), xlabel('w (rad)'), xlim([0 2*pi]), ylim([-0.2 0.2])
All of this becomes much simpler if the time domain signal, x[n], is causal, i.e., 0 for n < 0.
0 Kommentare
ZHU CHENG-CHIH
am 18 Apr. 2020
In Digital Signal Processing(DSP) class, we implemet it use matrix.
function [X] = dtft(x,n,w)
% Computes Discrete-time Fourier Transform
% [X] = dtft(x,n,w)
% X = DTFT values computed at w frequencies
% x = finite duration sequence over n
% n = sample position vector
% w = frequency location vector
X = exp(-1i*w'*n) * x.';
% X = x*exp(-j*n'*w);
end
Maybe you can get a DSP textbook to read it more. :D
We use this book : Digital Signal Processing Using MATLAB by Vinay K. Ingle, John G. Proakis.
There is a practice, Problems 3.1, in Chapter 3.
2 Kommentare
Jon
am 28 Jun. 2019
Yes - you can use the MATLAB FFT (fast fourier transform) function to compute DFT's. Please see the MATLAB documentation for detail https://www.mathworks.com/help/signal/ug/discrete-fourier-transform.html
1 Kommentar
Matt J
am 28 Jun. 2019
You could try using symsum in the Symbolic Math Toolbox. Why do you need a continuous-frequency result?
6 Kommentare
Matt J
am 29 Jun. 2019
Bearbeitet: Matt J
am 29 Jun. 2019
The DTFT is a sum of complex exponentials. My suggestion was that you might be able to compute that symbolically with symsum.
But before you do that, you should be sure you really need a symbolic, continuous-frequency result, instead of the discrete-frequency result that FFT already offers you.
SAMIR PAUL
am 13 Apr. 2023
Bearbeitet: SAMIR PAUL
am 13 Apr. 2023
clc;
b=[1];
a=[1,-0.5];
w=-pi:pi/100:pi;
[h]=freqz(b,a,w);
%magnitude spetrum
subplot(2,2,1);
plot(w/pi,abs(h));
xlabel('normalized frequency \omega/\pi');
ylabel('magnitude');
title('magnitude spetrum of the transfer function');
%phase spetrum
subplot(2,2,2);
plot(w/pi,angle(h));
xlabel('normalized frequency \omega/\pi');
ylabel('Phase in radians');
title('phase of the transfer function');
%real part
subplot(2,2,3);
plot(w/pi,real(h));
xlabel('normalized frequency \omega/\pi');
ylabel('magnitude');
title('real part of the transfer function');
%imaginary part
subplot(2,2,4);
plot(w/pi,imag(h));
xlabel('normalized frequency \omega/\pi');
ylabel('magnitude');
title('imaginary part of the transfer function');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1354138/image.jpeg)
0 Kommentare
Kavita Guddad
am 4 Jun. 2023
Bearbeitet: Kavita Guddad
am 6 Jun. 2023
clc;clear all;close all;
x=[2 3 6 -3];
n=0:length(x)-1;
w=-pi:0.01:pi;
[X] = dtft(x,w);
subplot(311);stem(x);title('Signal'); xlabel('time index');ylabel('amplitude');
subplot(312);plot(w,abs(X));title('Mag. plot'); xlabel('Frequency w in rad');ylabel('Mag.');
subplot(313);plot(w,angle(X));title('Phase plot'); xlabel('Frequency w in rad');ylabel('Phase angle');
function [X] = dtft(x,w)
% Computes Discrete-time Fourier Transform
% X = DTFT values computed at w frequencies
% x = finite duration sequence over n
% n = sample position vector
% w = frequency location vector
for i=0:length(w)-1
X(i+1)=0;
for k=0:length(x)-1
X(i+1) =X(i+1)+ exp(-1i*w(i+1)*k)* x(k+1);
end
end
end
%Note:you can change x and w vector values as per requirement
%e.g w=w=-3*pi:0.01:3*pi;;
2 Kommentare
Paul
am 4 Jun. 2023
On the first iteration through the for loop, i = 0 and X(i) = 0 will cause an error because array indices have to be positive integers or logical values. Also dtft is missing a closing "end" statement.
Kavita Guddad
am 6 Jun. 2023
Here is the corrected code.
Even n is not required
%Main program
clc; clear all; close all;
x=[2 3 6 -3];
n=0:length(x)-1;
w=-pi:0.01:pi;
X = dtft(x,w);
subplot(311);stem(n,x);title('Signal1'); xlabel('time index');ylabel('amplitude');
subplot(312);plot(w, abs(X));title('Mag. Plot');xlabel('Frequency w in radt');ylabel('Mag.');
subplot(313);plot(w, angle(X));title('Phase plot'); xlabel('Frequency w in rad');ylabel('amplitude');
%Function Program
function [X] = dtft(x,w)
% Computes Discrete-time Fourier Transform
% X = DTFT values computed at w frequencies
% x = finite duration sequence over n
% n = sample position vector
% w = frequency location vector
for i=0:length(w)-1
X(i+1)=0;
for k=0:length(x)-1
X(i+1) =X(i+1)+ exp(-1i*w(i+1)*k)* x(k+1);
end
end
end
Siehe auch
Kategorien
Mehr zu Transforms finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!