DFT的计算、FFT的基础代码、FFT的横纵坐标问题(matlab)

最新推荐文章于 2024-05-29 17:27:53 发布

安星眠

最新推荐文章于 2024-05-29 17:27:53 发布

阅读量1.1w

点赞数 12

文章标签：信号处理 matlab

本文链接：https://blog.csdn.net/qq_44827133/article/details/123142243

版权

FFT的定义

FFT：快速傅里叶变换，是DFT的快速算法。

DFT（Discrete Fourier Transform）：离散傅里叶变换。在DTFT之后，将傅里叶变换的结果也进行离散化，就是DFT。

DTFT（Discrete-time Fourier Transform）：离散时间傅里叶变换。这里的“离散时间”指的是时域上式离散的，也就是计算机进行了采样。不过傅里叶变换后的结果依然是连续的。

也就是说：DTFT时域离散、频域连续；DFT时域离散、频域离散。

FFT和DFT中的N

DFT数学计算公式：
$\begin{array}{l} X(k) = \sum\limits_{n = 0}^{N - 1} {x(n){e^{ - j\frac{{2\pi kn}}{N}}}} \\ x(n) = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {X(k){e^{j\frac{{2\pi kn}}{N}}}} \end{array}$
根据公式在matlab上自定义dft.m

function[Xk]=dft(xn,N);
n=[0:1:N-1];
k=[0:1:N-1];
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
end

然后设置

xn=[ ];
y1=fft(xn，N)
y2=dft(xn,N)

这个比较简单，结论为MATLAB中对N点序列做fft的结果与N点dft的结果是完全一样的，没有对N点序列进行补0后再做操作。

2.FFT中的n
Y = fft(X,n) 返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同。
如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n。
如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n。

3.DFT中的N
例：对R4(n)作不同N点数的DFT，分别为4点，8点，16点。
在这里插入图片描述
用fft( )函数之间计算也可，DFT和FFT两者计算后的结果是完全相同的。

xn=[1 1 1 1];
Y1 = fft(xn,4)
Y2 = fft(xn,8)

j
在这里插入图片描述
注：
作出它们的频谱图X(f)会发现，由于采样率默认相同，故它们的fs相同，不同的N带来的是频率间隔的不同。故其实它们的连续谱形状应该相同，但是分辨率不同导致了它们的差异。

FFT基础代码

一须知：
1.我们通常只研究信号FFT变换后形成的离散周期信号的一个周期，周期大小为采样频率Fs，且信号关于Fs/2对称。在FFT后的横坐标中最大不会超过Fs。
2.

T=1/Fs;
df=Fs/N;
t=(0:N-1)*T;

df为频率分辨率，相同的采样点数N下，fs越低则频谱分辨率越高，但会带来高频混叠的影响。x(n)有N个点，X(k)对应的也有N个点，在频域中第n个点对应的频率为（n-1)*df。
3.由于对称性，我们通常研究单侧频谱。当采样点数为N（通常为2的次幂法）时，单侧频谱序列个数如下：
当N为奇数时为（N+1)/2,
当N为偶数时为 (N/2)+1。
4.信号FFT变换后，当k=0，即第一个点，也是频率为0的这个点代表着信号的直流分量。有的信号需要去除直流分量，有一种方法：将其傅里叶变化后的第一个点的幅值设为0。
5.在双侧频谱转换成单侧频谱时，根据权系数的统一性，相应频谱(k≠0)单边谱的幅度谱幅值大小，是双边谱中的两倍。

二常用FFT代码

Fs=1000; %采样频率
L=1500;  %采样点数
T=1/Fs; %时间间隔
df=Fs/L;  %频率分辨率
t=(0:L-1)*T; % 时间序列
X = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t));
plot(1000*t(1:50),X(1:50)) %*1000是时间的单位换算
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')
Y = fft(X); %Y的结果是L个复数
P2 = abs(Y/L);% 也有写成abs(Y)/L，数学上是等价的，P2为双侧频谱的幅值
P1 = P2(1:L/2+1);%因为是偶数，所以单侧频谱包含点的个数为（L/2）+1
P1(2:end-1) = 2*P1(2:end-1);%除去k=0、k=L/2两点，其余点的幅值*2，P1的幅值就是对应频率真实的幅值。
f = (0:(L/2))*(Fs/L);
figure
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

稍作解释：
1.上一篇讲到，傅里叶的目的是将一个复杂信号分解成几个正弦信号之和，并在频域图上对分解后的正弦信号的两个指标（频率和幅度）进行展现。所以以上述信号为例，当你看到这个信号X = 0.7sin(2pi50t) + sin(2pi120t));你就应该知道它FFT变换后的图像，如图1。图2为未对FFT后纵坐标处理时的照片，可以看出5252/L=525*2/1500=0.7,正好差了2/L倍。
在这里插入图片描述

FFT的横纵坐标

1.先说纵坐标为什么*2/L，这个网上的说法太多了，我找了我认为两个靠谱的说法放在了后面参考文章中。值得注意的是并不是所有的点FFT后都需要乘，以L等于偶数同时离散时间序列为实数为例，第一个和第(L/2)+1个，也就是k=0、k=L/2时，需要×1/L。以L=8为例。如图L为偶数，所以单侧频谱的个数为(L/2)+1=5,第一个和第五个×1/L,其余×2/L,(途中P1已经是×1/L的结果了)。在这里插入图片描述
2.关于横坐标的理解见下图。