数字信号处理实验三 FFT 应用及 CZT （fft在快速卷积，相关，功率谱及CZT应用）

快速傅里叶变化FFT的应用

前言
快速傅里叶变换
快速卷积计算
快速相关计算
功率谱计算
线性调频Z变换（CZT）
全部程序可点此处下载

前言

傅里叶变换在时频域转换和频域分析上有着重要的作用。但是如果直接用DFT公式进行计算当序列长度N很大时，运算量大，且运算耗时。而快速傅里叶变换（FFT）就可以很好的解决这个问题。本次实验便是利用FFT算法，在快速卷积，快速自相关等理论基础上，使用matlab编写程序实现计算卷积，信号自相关和功率谱，实现CTZ变换。

快速傅里叶变换

N 点的 DFT 可以分解为两个 N/2 点的 DFT，每个 N/2 点的 DFT 又可以分解为两个 N/4 点的 DFT。依此类推，当 N 为 2 的整数次幂时由于每分解一次降低一阶幂次，所以通过 M 次的分解，最后全部成为 N/2 个 2 点 DFT 运算。以上就是按时间抽取的快速傅立叶变换(FFT)算法。当需要进行变换的序列的长度不是 2 的整数次方的时候，为了使用以 2 为基的 FFT，可以用末尾补零的方法，使其长度延长至 2 的整数次方。
在matlab中，fft()函数可以帮助我们方便的使用快速傅里叶变换，ifft()函数可以让我们方便的使用快速傅里叶逆变换。

快速卷积计算

快速卷积就是用快速傅里叶变换的方法来进行卷积。其方法与DFT实现卷积相同，只是将DFT改为FFT。应用时域频域卷积定理，时域卷积等于频域相乘。即两序列卷积等于两序列的FFT相乘后再将结果做IFFT即可得到时域卷积结果。过程图如下图所示。

例1
序列一：x(n)=cos(w n2 ) , 0<=n<=40 , w=0.025pi
序列二：h(n)=
[-0.9754, -0.7820, -0.4142, 0.092, 0.6871, 1.31287, 1.908,
2.4142, 2.7820, 2.9754, 2.9754, 2.7820, 2.4142, 1.908,
1.31287, 0.6871, 0.092, -0.4142, -0.7820, -0.9754]
对这两个序列使用快速傅里叶的方法进行卷积运算。在matlab中编写程序，并使用直接卷积的方法进行验证。
使用freqz函数求频响特性，再使用abs和angle分别求幅频响应和相频响应。
程序

clear;close all;
n=0:1:40;
w=0.025*pi;
x=cos(w*(n.^2));
h=[-0.9754,-0.7820,-0.4142,0.092,0.6871,1.31287,1.908,2.4142,2.7820,2.9754,2.9754,2.7820,2.4142,1.908,1.31287,0.6871,0.092,-0.4142,-0.7820,-0.9754];
N1=length(x);
N2=length(h);
N=N1+N2-1;        %计算卷积序列长度
X=fft(x,N);     %计算两个序列fft
H=fft(h,N);
Y=X.*H;      %频域相乘
y=ifft(Y);
m=0:1:N-1;
subplot(211);
stem(m,y);
title('fft法');
y1=conv(x,h);  %卷积函数验证
subplot(212)
stem(m,y1);
title('直接卷积法');

运行结果如下

快速相关计算

快速相关即使用FFT的方法计算相关函数。实际上应用的是相关定理，计算流程如下图所示。

计算自相关时，两个序列相同，计算出rxx就为自相关序列，Rxx为其功率谱。

例
序列x(n)=[-0.9754, -0.7820, -0.4142, 0.092, 0.6871, 1.31287, 1.908,
2.4142, 2.7820, 2.9754, 2.9754, 2.7820, 2.4142, 1.908,
1.31287, 0.6871, 0.092, -0.4142, -0.7820, -0.9754]
要求计算 h(n)的自相关函数，并将结果显示出来。
程序如下

clear;close all;
x=[-0.9754,-0.7820,-0.4142,0.092,0.6871,1.31287,1.908,2.4142,2.7820,2.9754,2.9754,2.7820,2.4142,1.908,1.31287,0.6871,0.092,-0.4142,-0.7820,-0.9754];
N=length(x); %fft法
xn=[x,zeros(1,N-1)];  %补零
X=fft(xn);           %正变换
R=X.*conj(X);       %套用公式
r=ifft(R);            %逆变换
r1= circshift(r, N-1); %移位r0=xcorr(x,x);     %自带函数法
n=1-N:N-1;
subplot(211);
plot(n,r1);
subplot(212);
plot(n,r0);

分别使用fft法和使用matlab自带相关函数进行计算，计算结果相同，证明fft方法正确性。运行结果如下

功率谱计算

例
信号设为：x(t)=cos2πf1t + cos2πf2t + cos2πf3t
f1=42Hz, f2=45Hz, f3=47Hz,采样频率 fs=200Hz,采样点数 N=128。
计算信号的功率谱，将功率谱图显示出来。
计算方法与计算自相关方法类似，但自相关序列的傅里叶变换就是序列的功率谱，因此，计算出Rxx即可。
程序如下

close all;
clear;
%抽样
f1=42;f2=45;f3=47;
N=128;T=1/200;
t=0:T:(N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t);
%计算自相关序列，方法和第二步相同
xn=[x,zeros(1,N-1)];
X=fft(xn);
R=X.*conj(X);%自相关函数频谱即为功率谱
plot(abs(R))

运行结果如下

线性调频Z变换（CZT）

CZT 是 Z 平面上的一段螺线等间隔采样，，可以利用 FFT 计算 CZT。
其计算公式如下。

或使用如下流程图，利用fft计算
例
信号设为：x(t)=cos2πf1t + cos2πf2t + cos2πf3t
f1=42Hz, f2=45Hz, f3=47Hz
要求
a. 用 FFT 计算信号的频谱，并显示信号的幅频和相频波形，坐标横轴以频率 Hz为单位。
b. 以下列参数计算 CZT：
起始频率：40Hz，点数：20，频率间隔：0.5Hz
将 CZT 的幅频、相频图显示在屏幕上，坐标横轴以频率 Hz 为单位。

在matlab中，可以使用函数czt(x,M,w,a)来计算CZT变换。利用如下公式计算出w和a，M为取样点数，带入函数中，即可球的信号的CZT

程序如下

close all;
clear;
%抽样
f1=42;f2=45;f3=47;
N=128;fs=200;Ts=1/fs;t=(0:N-1)*Ts;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t);
X=fft(x);     %快速傅里叶变换
f=(0:N-1)*fs/N;  %确定频率范围
subplot(411);
plot(f,abs(X));
xlabel('频率HZ');
title('fft幅度');
subplot(412);
plot(f,angle(X));
xlabel('频率HZ');
title('fft相位');fs1=0.5;
fl=40;
M=20;
fh=fl+M*fs1;    %确定抽取频率范围
w=exp(-1j*2*pi*(fh-fl)/(fs*M)); %计算参数Wn
a=exp(1j*2*pi*fl/fs);  %计算参数An
xk=czt(x,M,w,a);   %代入计算czt结果m=0:1:M-1;
f0=fs1*m+fl;    %确定频率范围
subplot(413);
plot(f0,abs(xk));
title('czt幅度');
xlabel('频率HZ');
subplot(414);
plot(f0,angle(xk));
title('czt相位');
xlabel('频率HZ');

运行结果如下

全部程序可点此处下载