问题引入
上面是DFT和IDFT的公式，IDFT先不谈。在实践中常使用FFT算法来快速算出DFT，获得时域采样信号x(n)的频率幅度谱。
例如：
1
2
| N=1500;
xn=0.5*sin(2*pi*75/N*n)+3*cos(2*pi*40/N*n)+0.7*sin(2*pi*350/N*n);
|
上面这个时域采样信号xn由三个频率叠加而成，我们知道频率的概念是1秒信号经过了多少个周期，
那么我们假设采样率是fs=1000，共采样N=1500点，那么0.5*sin(2*pi*75/N*n)的频率即为75*fs/N=50，因此第一个信号频率为50。
通过Matlab进行FFT，可以画出这样一张频率幅度谱：
这里就引申出我们的标题：DFT为什么能求频率幅度谱？DFT后的X[k]与x(n)幅度的关系？DFT/IDFT底层数学原理？
对于这个问题，我看了不少网上的博客、视频、书籍，总觉得讲的不清不楚，有人说X[k]的结果就是x(n)的幅度，有人说X[k]的结果除以N就是x(n)的幅度，有人说X[k]的结果除以N/2就是x(n)的幅度，就没有一个比较清楚的证明到底X[k]的值和x(n)幅度的关系。
如果你也对此有疑问，现在我证明给你看！
铺垫一些小公式
要看懂我的证明，首先需要铺垫一些小公式，防止你不知道！
DFT公式证明
我将公式的证明整理成了PPT，下面我直接复制，在证明中使用到了上面铺垫的小公式，注意关联起来看就可以理解了！
DFT公式分解为4部分
先考虑k1=0的情况:
再考虑k1≠0的情况:
DFT计算后，X(k)与x(n)的关系：
可以看到，DFT公式本质上可以完全拆成三角函数的计算，在此过程中利用
高中的积化和差公式和三角函数在若干个周期的累加为0，
可以得到X[k]和原信号x(n)幅度的关系。
Matlab FFT示例代码
下面我们给出使用FFT对x(n)进行DFT计算后，绘制频率幅度谱的Matlab代码，看官可以发现，这上述公式的X(k)结果是完全对应的！
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
| %功能说明：
%生成时域采样信号0.5sin(2pi*k/N*n)
%通过fft转化到频域，绘制频率幅度谱，横坐标是频率f=k*fs/N，纵坐标是实际幅度
%再通过ifft转化回时域，绘制ifft结果信号，和原信号进行对比
%采样率1000
fs=1000;
%在N个采样点中，一共振动了k个周期
k=75;
%采样总点数
N=1500;
%n表示采样的第n个点
n=0:N-1;
%xn是采样得到的x(n)时域信号
xn=0.5*sin(2*pi*k/N*n)+3*cos(2*pi*40/N*n)+0.7*sin(2*pi*350/N*n);
xk=fft(xn);
%abs是求模
P2=abs(xk/N);
P1=P2(1:N/2+1);
P1(2:end-1)=2*P1(2:end-1);
%半频谱 f=k*fs/N
f=fs*(0:(N/2))/N;
plot(f,P1)
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")
|
IDFT公式证明
Matlab调用FFT/IFFT并绘图
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
| %功能说明：
%生成时域采样信号0.5sin(2pi*k/N*n)
%通过fft转化到频域，绘制频率幅度谱，横坐标是频率f=k*fs/N，纵坐标是实际幅度
%再通过ifft转化回时域，绘制ifft结果信号，和原信号进行对比
%采样率1000
fs=1000;
%在N个采样点中，一共振动了k个周期
k=75;
%采样总点数
N=1500;
%n表示采样的第n个点
n=0:N-1;
%xn是采样得到的x(n)时域信号
xn=0.5*sin(2*pi*k/N*n)+3*cos(2*pi*40/N*n)+0.7*sin(2*pi*350/N*n);
xk=fft(xn);
%abs是求模
P2=abs(xk/N);
P1=P2(1:N/2+1);
P1(2:end-1)=2*P1(2:end-1);
%半频谱 f=k*fs/N
f=fs*(0:(N/2))/N;
plot(f,P1)
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")
% 反傅里叶变换
X = ifft(xk, 'symmetric');
% 绘制输入信号和还原信号的图像
figure
subplot(2,1,1)
plot(n, xn)
title('Input Signal')
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,1,2)
plot(n, X)
title('Signal after IFFT')
xlabel('Time (s)')
ylabel('Amplitude')
|
