【 MATLAB 】使用 MATLAB 作圖討論有限長序列的 N 點 DFT（強烈推薦）（含MATLAB指令碼）
阿新 • • 發佈：2018-12-18
但是這篇博文我最初設計的是使用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 ');
當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上的取樣間隔，也就是頻率解析度。
補零給出了一種高密度的譜並對畫圖提供了一種更好的展現形式。
但是它並沒有給出一個高解析度的譜，因為沒有任何新的資訊附加到這個訊號上；而僅是在資料中添加了額外的零值。
下篇博文我們繼續說明，如何才能得到高解析度的譜，我們的方法是從實驗中或觀察中獲得更多的資料。