文章目录
- 问题引入
- 铺垫一些小公式
- DFT公式证明
- DFT公式分解为4部分
- 先考虑k1=0的情况:
- 再考虑k1≠0的情况:
- DFT计算后,X(k)与x(n)的关系:
- Matlab FFT示例代码
- IDFT公式证明
- Matlab调用FFT/IFFT并绘图
问题引入
上面是DFT和IDFT的公式,IDFT先不谈。在实践中常使用FFT算法来快速算出DFT,获得时域采样信号x(n)的频率幅度谱。
例如:
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)结果是完全对应的!
%功能说明:
%生成时域采样信号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并绘图
%功能说明:
%生成时域采样信号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')