信号处理——信号频域变换

作者：桂。

时间：2017-03-07  18:40:05

链接：http://www.cnblogs.com/xingshansi/p/6516256.html 

声明：欢迎转载，不过记得注明出处哦~

前言

 本文主要介绍信号频域变换的常用方式，主要介绍基本的三种方法，实现时域数字信号的傅里叶变换：

　　1）基本DFT变换；

　　2）FFT变换

　　3）逆序级联FFT

三种方法实现方式略有差别，但本质完全相同。为了方便大家理解，本文给出基本MATLAB实现。内容的理论部分，多有借鉴他人，相应链接在最后一并给出。

 〇、写在前面

关于信号时域连续、离散以及各自频域的对应关系，可以参考之前的一篇博文。为了便于理解，此处给一个DTFT—>DFT的解释。

MATLAB代码：

%DTFT
clc;clear all;close all
N=8;                          %原离散信号有8点
n=[0:1:N-1]    ;              %原信号是1行8列的矩阵
xn=0.5.^n;                    %构建原始信号，为指数信号

w0=[-2*N:.1:2*N]*2*pi/N;       %方便观察，此处取4个周期-4*pi~4*pi  
X0=xn*exp(-j*(n'*w0));         %求dtft变换，采用原始定义的方法，对复指数分量求和而得
subplot(221)
plot(n,xn,'r');hold on;xlabel('时间')
stem(n,xn,'k');grid on;
title('时域采样后的信号(指数信号)');ylabel('幅值')
subplot(222);
plot(w0,abs(X0),'k');ylabel('幅值');xlabel('角频率');%逼近连续，其实这么处理不正确，仅仅是方便理解
title('DTFT变换');
xlim([-4*pi,4*pi]);grid on;

%DFT
w=[-2*N:1:2*N]*2*pi/8;      
X=xn*exp(-j*(n'*w));         %对DTFT进行频域采样
subplot(223)
stem(n,xn,'k');ylabel('幅值');xlabel('时间');grid on;
title('频域采样后的信号(指数信号)');
subplot(224);
plot(w0,abs(X0),'r');hold on;
stem(w,abs(X),'k');ylabel('幅值');xlabel('角频率')
title('DFT变换');
xlim([-4*pi,4*pi]);grid on;

　　对应结果图：

一、基本DFT变换

DFT变换有：

$F(k) = \sum^{N-1}_{n=0} f(n) e^{\frac{-j2\pi kn}{N}}$

写成矩阵的形式：

对应的转换为VanderMonde矩阵，对应的MATLAB代码：

clc;clear all;close all;
fs = 200;
f0 = 30;
t = 0:1/fs:2;
x = sin(2*pi*f0*t');
N = length(x);
%DFT
X_DFT = exp(-j*2*pi/N).^([0:N-1]'*[0:N-1])*x;
%绘图
xf = linspace(0,fs,N);
plot(xf,abs(X_DFT),'k');grid on;
xlabel('频率（Hz）');
ylabel('幅度');
title('DFT结果');

　　也可以将DFT的代码换为Vander矩阵形式，只需修改一句：

% X_DFT1 = exp(-j*2*pi/N).^([0:N-1]'*[0:N-1])*x;
%修改为下句
X_DFT = fliplr(vander(exp(-j*2*pi/N).^([0:N-1]))).'*x;

结果图：

二、FFT变换

 DFT到FFT常用的有奇偶抽取/前后抽取，两种实现方式，原理此处不再堆砌，可以参考维基百科：

给出对应的MATLAB代码：

clc;clear all;close all;
fs = 200;
f0 = 30;
t = 0:1/fs:2;
x = sin(2*pi*f0*t');
N = length(x);
%FFT
X_FFT = fft(x);
%绘图
xf = linspace(0,fs,N);
plot(xf,abs(X_FFT),'k');grid on;
xlabel('频率（Hz）');
ylabel('幅度');
title('FFT结果');

结果图：

如果严格按照FFT定义，则应该将序列补为$2^N$长度（事实上，内置fft自动完成此操作），代码修改为：

clc;clear all;close all;
fs = 200;
f0 = 30;
t = 0:1/fs:2;
x = sin(2*pi*f0*t');
N = length(x);
%FFT
Len = 2^nextpow2(length(x));
X_FFT = fft(x,Len);
%绘图
xf = linspace(0,fs,Len);
plot(xf,abs(X_FFT),'k');grid on;
xlabel('频率（Hz）');
ylabel('幅度');
title('FFT结果');

　　对应结果图：

三、逆序级联FFT

为了方便，假设信号$f(n)$长度为2的整次幂。对于数据量过大的情况，可以通过串/并的不断转化，对数据进行分流，从而利用多个芯片实现一个大芯片的任务（此处为个人理解，仅供参考）。给出对应的流图：

以及具体的推导公式：

步骤可以简化为：分段FFT，即公式[.]部分—>因子补偿，即公式{.}部分—>分段FFT，即公式最外层部分。

给出具体的实现代码：

clc;clear all;close all;
fs = 200;
f0 = 30;
t = 0:1/fs:2;
x = sin(2*pi*f0*t');
Na = 2^nextpow2(length(x));
sig = [x;zeros(Na-length(x),1)];
N=16;M=Na/N;
sr=zeros(M,N);
for m=1:M
    sr(m,:)=sig((m-1)*N+1:m*N,1);
end
for n=1:N
    sr(:,n)=fft(sr(:,n));
end
P=zeros(M,N);
for p=0:M-1
    for n=0:N-1
        P(p+1,n+1)=exp(-j*2*pi*n*p/M/N);
    end
end
sr=sr.*P;
for m=1:M
    sr(m,:)=fft(sr(m,:));
end
%save result
X_CaFFT = zeros(length(sig),1);
for n=1:N
    X_CaFFT((n-1)*M+1:n*M,1)=sr(:,n);
end
%绘图
xf = linspace(0,fs,Na);
plot(xf,abs(X_CaFFT),'k');grid on;hold on;
xlabel('频率（Hz）');
ylabel('幅度');
title('CaFFT结果');

　　对应结果图：

四、三种算法对比分析

三种结果效果完全等价，分析其各自性能。

对于快速运算中，一个蝶形运算，对应一个乘法，两次加法。

  DFT FFT CaFFT
加法器 $N^2$ $N log_2N$ $N log_2N+N$
乘法器 $N(N-1)$ $2Nlog_2N$ $2Nlog_2N$

可见FFT所需资源最少，但如果数据量过大导致直接FFT，硬件无法满足条件，则CaFFT是最优的选择。

 

参考：

郑君里：《信号与系统》第二版;

张大炜：一种新的级联FFT算法.

posted @ 2017-03-07 22:18  LeeLIn。  阅读(5380)  评论(0)  编辑  收藏  举报