DFT、DTFT和DFS是否已經把你搞暈

http://m.sohu.com/a/249761071_119709 缺失的圖5和6在 http://www.cnblogs.com/BitArt/archive/2012/11/24/2786390.html

一幅圖弄清DFT與DTFT,DFS的關係http://www.cnblogs.com/iamhcc/p/5967884.html

https://blog.csdn.net/u012676933/article/details/53141117 假設採樣頻率爲Fs，採樣點數爲N，作FFT以後，某一點n（n從1開始）表示的頻率爲：Fn=(n-1)*Fs/N；該點的模值除以N/2就是對應

%該頻率下的信號的幅度（對於直流信號是除以N），具體網頁內容見下：

很基礎的問題每每很重要，作仿真時候有一個點的差錯都會致使結果的錯誤。在網上找了前人寫的東東，總結下但願對你們有幫助，讓你們少走一些彎路。

1.  信號的時域採樣點N和頻域採樣點數相同

%##################################################################

clear all; close all;
Adc = 1.25;                   %直流份量幅度
A1 = 1;                       %頻率F1信號的幅度
A2 = 0.25;                     %頻率F2信號的幅度
F1 = 100;                    %信號1頻率(Hz)
F2 = 1000;                    %信號2頻率(Hz)
Fs = 5120;                   %採樣頻率(Hz)

Ts = 1/Fs;                     %時域採樣間隔，採樣週期
P1 = -30;                    %信號1相位(度)
P2 = 90;                     %信號2相位(度)
N = 256;                      %採樣點數，請注意時域採樣N點，FFT時也是N點！
t = [0 : 1/Fs : N/Fs];        %採樣時刻，注意不是序列是真正的時間！t = [0:N]*Ts

% t = [0 : N]//Fs;  %共採樣257個點，爲了好看起見多采樣了最後一個點，最後一個點在實際中應該是下一個採樣週期的第一個點

%生成信號
%S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180); %考慮相位
signal1=A1*sin(2*pi*F1*t); % 週期T1 = 2*pi/(2*pi*F1) = 0.01s， 採樣週期Ts = 0.00019531 < (T1)/2

signal2=A2*sin(2*pi*F2*t); % 週期T2 = 2*pi/(2*pi*F2) = 0.001s，採樣週期Ts = 0.00019531 < (T2)/2
S=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t);                       %無直流偏置
S=Adc+A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t);                   %直流偏置
%S=Adc+[A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t)]/2;

S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180); %考慮相位
%顯示原始信號
plot(S);
title('原始信號');

figure;
Y = fft(S,N);            %作FFT變換
Ayy = (abs(Y));          %取模
plot(Ayy(1:N));          %顯示原始的FFT模值結果
title('FFT 模值');

F = ([1:N] - 1) * Fs / N;               %換算成實際的頻率值， N*Fs/2 對應着w = pi
% 由於MATLAB中FFT的變換矩陣不是一個酉矩陣(Unitary Matrix)，該陣除以1/sqrt(N)就是個酉矩陣。故通過變換後對信號有放大做用，

% 因此要在fft處理後結果除以N/2來，修正此「放大」做用。可是結果在直流的那一點是錯誤的，實際上直流應該除以N修正。

Ytemp = 2/N * fft(S,N); % 作FFT變換後除以N/2修正

figure; plot(F, abs(Ytemp));title('實際幅度-頻率曲線圖 --- abs前修正');

figure;
Ayy = Ayy / (N / 2);                    %換算成實際的幅度
Ayy(1) = Ayy(1) / 2;                    %能看出來爲何直流信號除以N了嗎？我開始沒看出來啊

%假設採樣頻率爲Fs，採樣點數爲N，作FFT以後，某一點n（n從1開始）表示的頻率爲：Fn=(n-1)*Fs/N；該點的模值除以N/2就是對應

%該頻率下的信號的幅度（對於直流信號是除以N）

 

根據公式驗證 --- 某點n所表示的頻率爲： Fn=(n-1)*Fs/N，由於沒有fftshift

                                     f1=（6-1）*5120/256=100 （Hz） ---- 驗證結果正確

                                     f2=（51-1）*5120/256=1000（Hz）---- 驗證結果正確

 根據公式驗證 --- 對於n=1點的信號，是直流份量，幅度即爲A1/N   --- A1=320/256=1.25 --- 驗證結果正確

                 對於n點（n≠1，且n<=N/2）  幅度A --- A=An/(N/2)=128/（256/2）=1 --- 驗證結果正確

plot(F(1:N/2), Ayy(1:N/2));             %顯示換算後的FFT模值結果， 只看0~N*Fs/2正頻譜部分
title('實際幅度-頻率曲線圖 --- abs後修正');

figure;

% 相位的計算可用函數atan2(b,a)計算。atan2(b,a)是求座標爲(a,b)點的角度值，範圍從-pi到pi。

% atan2(500, 148)=x,結果是弧度，換算爲角度就是180*(-x)/pi=相位值。

Pyy = [1 : N/2];
for i = 1 : N/2
Pyy(i) = phase(Y(i));                %計算相位
Pyy(i) = Pyy(i) * 180 / pi;          %換算爲角度
end;
plot(F(1 : N/2), Pyy(1 : N/2));      %顯示相位圖
title('相位-頻率曲線圖');

 

根據FFT結果以及上面的分析計算，--- 頻率,幅度.相位....咱們就能夠寫出信號的表達式了，它就是咱們開始提供的信號。------ 信號重建

2.  信號的時域採樣點N和頻域採樣點數NFFT不一樣

%若是採樣的時間週期是信號週期的倍數，可能泄露就會避免。頻譜泄露是因爲非整週期採樣致使的。

% 我之前有個錯誤的觀念：信號的時域採樣點N必須和FFT的計算點數NFFT相同，纔會給處理和解釋帶來便利。

% 原來是模型的不一樣產生的影響，採用模擬頻率 f 建模不會產生末尾補零使得FFT頻率不一致的問題，看下面的程序。

clear; close all;

% f1 = 20Hz, f2 = 40Hz %若 f1 = 15 能夠觀察頻譜泄露狀況

%  （1）數據個數Ndata=32，FFT所用的採樣點數NFFT=32
fs=240; %採樣頻率
Ndata=32; %數據長度
N=32; %FFT的數據長度
n=0:Ndata-1;t=n/fs;   %數據對應的時間序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);   %時間域信號， x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); 能夠發現頻譜泄露
y=fft(x,N);   %信號的Fourier變換
mag=abs(y);    %求取振幅
f=(0:N-1)*fs/N; %真實頻率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率以前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');grid on;

% （2）數據個數Ndata=32，FFT所用的採樣點數NFFT=126
Ndata=32;   %數據個數
N=128;     %FFT採用的數據長度
n=0:Ndata-1;t=n/fs;   %時間序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真實頻率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率以前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');grid on;

% （3）數據個數Ndata = 136，FFT所用的採樣點數NFFT=128
Ndata=136;   %數據個數
N=128;     %FFT採用的數據個數
n=0:Ndata-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真實頻率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率以前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');grid on;

% （4）數據個數Ndata=136，FFT所用的採樣點數NFFT=512
Ndata=136;    %數據個數
N=512;    %FFT所用的數據個數
n=0:Ndata-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真實頻率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率以前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
% （5） 數據個數Ndata=136，FFT所用的採樣點數NFFT=4096
Ndata=136;    %數據個數
N=4096;    %FFT所用的數據個數
n=0:Ndata-1;t=n/fs; %時間序列
x=0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真實頻率

figure,plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率以前的振幅
xlabel('頻率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;

3.  信號的時域採樣點N和頻域採樣點數NFFT不一樣 -- 比較兩種不一樣的信號模型

clear; close all;

% 該信號的數字週期 N = 8，模擬週期T = N*Ts = 0.008s ，實際頻率 f = 125 Hz

w = pi/4;
N = 2*pi/w; % N = 8 
n = 0:N-1; % n = 0:N 可能會更好看一些，可是要清楚第N+1點但是下一個採樣週期的第一個點
x = sin(n*w);  % 不含有采樣信息的x， 補零的FFT對頻譜有影響
h=plot(n, x, '-o'); %注意n沒有定標，沒有物理含義！
Fs = 1000; %採樣頻率爲1000Hz
Ts = 1/Fs; 
t = n*Ts; % t 時間序列的給法永遠只有這麼一種，請銘記！
T = N*Ts; %模擬週期T
f = w*Fs/(2*pi); % 信號的真實頻率f
figure; plot(t, x)
% freq = n*Fs/N - Fs/2; %頻率序列的定標，注意此式子是原本面目，不要給此式子穿個馬甲就不認識了，下面是此式子的變形
% freq = (n/N - 1/2)*Fs;        freq = (n - N/2)/N * Fs; % 注意0~N-1時 第N 點對應着 Fs， N/2 點對應着 Fs/2， Fs/N 就是頻率分辨率
% freq = (n - N/2)/N;  就是歸一化頻率，從-0.5 到 0.5 - 1/N , 看看不是 -0.5 ~ 0.5 哦，實際結合硬件編程就要這樣幹

Nfft = 1024; % 此處 Nffft/N = 128，因此要對結果的頻譜進行128點的抽取，才能夠獲得正確的頻率定標！
X_Nfft = fftshift(abs(fft(x,Nfft)));
freqNormalized = ((0:Nfft-1) - Nfft/2)/Nfft;
freq_ = freqNormalized*Fs;
freq = freqNormalized(1:128:end)*Fs; % 結果和上面的結果徹底相同 -500  -375  -250  -125     0   125   250   375

X_ = X_Nfft(1:128:end);

figure; plot(freq_, X_Nfft);grid

figure; plot(freq , X_);grid

% -------------------------------------------------------------------------
Ndata=136;    %數據個數
N=4096;    %FFT所用的數據個數
n=0:Ndata-1;t=n/fs; %時間序列
x = 0.5*sin(2*pi*20*t)+2*sin(2*pi*40*t); %含有采樣信息， 補零後FFT對信號的頻譜沒影響
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N;   %真實頻率

figure,plot(f(1:N/2),mag(1:N/2)*2/N); %繪出Nyquist頻率以前的振幅 xlabel('頻率/Hz');ylabel('振幅'); title('Ndata=136 Nfft=512');grid on;