【 MATLAB 】使用 MATLAB 作图讨论有限长序列的 N 点 DFT(强烈推荐)(含MATLAB脚本)
这篇博文本来是和上篇博文一起写的:【 MATLAB 】离散傅里叶级数(DFS)与DFT、DTFT及 z变换之间的关系
但是这篇博文我最初设计的是使用MATLAB脚本和图像来讨论的,而上篇博文全是公式,因此,还是单独成立了一篇,但是我还是希望看这篇博文之前还是先看看上篇博文。
我默认你已经看了上篇博文。
本博文的讨论建立在一个案例的基础上:
设x(n)是4点序列为:
计算x(n)的4点DFT(N =4),以及(N = 8,N = 16, N = 128)
我们知道DFS是在DTFT上的频域采样,采样间隔为 2*pi / N.
DFT 是一个周期的DFS,因此DFT也是在DTFT上的频域采样,采样间隔同样为 2 * pi / N.
x(n)的波形图为:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
stem(n,x);
-
title('Discrete sequence');
-
xlabel('n');ylabel('x(n)');
-
xlim([0,5]);ylim([-0.2,1.2]);
我们先作出 x(n) 的 DTFT,也即是离散时间傅里叶变换,给出脚本和图像:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
stem(n,x);
-
title('Discrete sequence');
-
xlabel('n');ylabel('x(n)');
-
xlim([0,5]);ylim([-0.2,1.2]);
-
-
%Discrete-time Fourier Transform
-
K = 500;
-
k = 0:1:K;
-
w = 2*pi*k/K; %plot DTFT in [0,2pi];
-
X = x*exp(-j*n'*w);
-
% w = [-fliplr(w),w(2:K+1)]; %plot DTFT in [-pi,pi]
-
% X = [fliplr(X),X(2:K+1)]; %plot DTFT in [-pi,pi]
-
-
magX = abs(X);
-
angX = angle(X)*180/pi;
-
-
figure
-
subplot(2,1,1);
-
plot(w/pi,magX);
-
title('Discrete-time Fourier Transform in Magnitude Part');
-
xlabel('w in pi units');ylabel('Magnitude of X');
-
subplot(2,1,2);
-
plot(w/pi,angX);
-
title('Discrete-time Fourier Transform in Phase Part');
-
xlabel('w in pi units');ylabel('Phase of X ');
-
关于DTFT的编程案例,之前博文也有:【 MATLAB 】模拟信号采样及离散时间傅里叶变换(DTFT)案例分析
当N = 4时候的DFT,为了显示出DFT和DTFT之间的关系,我们把DTFT图形和DFT画到了同一张图上:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
-
-
%Discrete-time Fourier Transform
-
K = 500;
-
k = 0:1:K;
-
w = 2*pi*k/K; %plot DTFT in [0,2pi];
-
X = x*exp(-j*n'*w);
-
-
magX = abs(X);
-
angX = angle(X)*180/pi;
-
-
-
-
%DFT in N = 4
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 4');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 4');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
可见,DFT是DTFT上的等间隔采样点。
当N = 8时候的DFT,为了显示出DFT和DTFT之间的关系,我们把DTFT图形和DFT画到了同一张图上:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
-
-
%Discrete-time Fourier Transform
-
K = 500;
-
k = 0:1:K;
-
w = 2*pi*k/K; %plot DTFT in [0,2pi];
-
X = x*exp(-j*n'*w);
-
-
magX = abs(X);
-
angX = angle(X)*180/pi;
-
-
-
-
% Zero padding into N = 8 and extended cycle
-
N = 8;
-
x = [1,1,1,1,zeros(1,4)];
-
-
%DFT in N = 8
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 8');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 8');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
相比于N =4,N =8采样点数更密集了。
当N = 16时候的DFT,为了显示出DFT和DTFT之间的关系,我们把DTFT图形和DFT画到了同一张图上:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
-
-
%Discrete-time Fourier Transform
-
K = 500;
-
k = 0:1:K;
-
w = 2*pi*k/K; %plot DTFT in [0,2pi];
-
X = x*exp(-j*n'*w);
-
-
magX = abs(X);
-
angX = angle(X)*180/pi;
-
-
-
-
% Zero padding into N = 16 and extended cycle
-
N = 16;
-
x = [1,1,1,1,zeros(1,12)];
-
-
%DFT in N = 16
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 16');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 16');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
N = 128的情况:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
-
-
%Discrete-time Fourier Transform
-
K = 500;
-
k = 0:1:K;
-
w = 2*pi*k/K; %plot DTFT in [0,2pi];
-
X = x*exp(-j*n'*w);
-
-
magX = abs(X);
-
angX = angle(X)*180/pi;
-
-
-
% Zero padding into N = 128 and extended cycle
-
N = 128;
-
x = [1,1,1,1,zeros(1,124)];
-
-
%DFT in N = 128
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 128');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 128');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
-
可见,随着周期N 的增大,采样点越密集。
上面的程序都是我在我写的一个大程序中抽取出来的,最后给出所有程序的一个集合,也就是我最初写的一个程序:
-
clc;clear;close all;
-
-
% x(n)
-
N = 4;
-
n = 0:N-1;
-
x = [1,1,1,1];
-
-
stem(n,x);
-
title('Discrete sequence');
-
xlabel('n');ylabel('x(n)');
-
xlim([0,5]);ylim([-0.2,1.2]);
-
-
%Discrete-time Fourier Transform
-
K = 500;
-
k = 0:1:K;
-
w = 2*pi*k/K; %plot DTFT in [0,2pi];
-
X = x*exp(-j*n'*w);
-
% w = [-fliplr(w),w(2:K+1)]; %plot DTFT in [-pi,pi]
-
% X = [fliplr(X),X(2:K+1)]; %plot DTFT in [-pi,pi]
-
-
magX = abs(X);
-
angX = angle(X)*180/pi;
-
-
figure
-
subplot(2,1,1);
-
plot(w/pi,magX);
-
title('Discrete-time Fourier Transform in Magnitude Part');
-
xlabel('w in pi units');ylabel('Magnitude of X');
-
subplot(2,1,2);
-
plot(w/pi,angX);
-
title('Discrete-time Fourier Transform in Phase Part');
-
xlabel('w in pi units');ylabel('Phase of X ');
-
-
-
%DFT in N = 4
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 4');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 4');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
-
-
% Zero padding into N = 8 and extended cycle
-
N = 8;
-
x = [1,1,1,1,zeros(1,4)];
-
-
%DFT in N = 8
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 8');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 8');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
-
-
-
% Zero padding into N = 16 and extended cycle
-
N = 16;
-
x = [1,1,1,1,zeros(1,12)];
-
-
%DFT in N = 16
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 16');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 16');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
-
-
-
% Zero padding into N = 128 and extended cycle
-
N = 128;
-
x = [1,1,1,1,zeros(1,124)];
-
-
%DFT in N = 128
-
k = 0:N-1;
-
Xk = dft(x,N);
-
magXk = abs(Xk);
-
phaXk = angle(Xk) * 180/pi; %angle 。
-
k = 2*pi*k/N; % k in DFT is equal to 2*pi*k/N in DTFT
-
-
figure
-
subplot(2,1,1);
-
stem(k/pi,magXk);
-
title('DFT Magnitude of x(n) when N = 128');
-
xlabel('k');ylabel('Magnitude Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,magX,'g');
-
hold off
-
-
subplot(2,1,2);
-
stem(k/pi,phaXk);
-
title('DFT Phase of x(n) when N = 128');
-
xlabel('k');ylabel('Phase Part');
-
-
% Auxiliary mapping, highlighting the relationship between DFT and DTFT
-
hold on
-
plot(w/pi,angX,'g');
-
hold off
-
-
最后需要说明的是上面通过补零的方式来增大最初给出的有限长序列的周期,这样只会更加DFT在DTFT上的采样间隔,也就是频率分辨率。
补零给出了一种高密度的谱并对画图提供了一种更好的展现形式。
但是它并没有给出一个高分辨率的谱,因为没有任何新的信息附加到这个信号上;而仅是在数据中添加了额外的零值。
下篇博文我们继续说明,如何才能得到高分辨率的谱,我们的方法是从实验中或观察中获得更多的数据。
文章来源: reborn.blog.csdn.net,作者:李锐博恩,版权归原作者所有,如需转载,请联系作者。
原文链接:reborn.blog.csdn.net/article/details/83476742
- 点赞
- 收藏
- 关注作者
评论(0)