MATLAB实现数字信号处理范文合集

第一篇：MATLAB实现数字信号处理

数字信号处理

说 明 书

目录

一．摘要…………………………………3 二．课程设计目的………………………3 三．设计内容……………………………3 四．设计原理……………………………4 4.1．语音信号的采集…………………………….4 4.2．滤波器……………………………………….4 4．21．IIR滤波器原理…………………………………….4 4.22．FIR滤波器原理………………………………………5 五．设计步骤……………………………6 5.1录制女音………………………………………6 5.2采样语音信号并画出时域波形和频谱图……7 5.3采用双线性变换法设计IIR滤波器…………10 5.4窗函数法设计FFR滤波器………………......12 5.5用IIR滤波器对信号进行滤波………………14 5.6用FIR滤波器对信号进行滤波………………16 5.7男女声语音信号频谱特点分析………………19 5.8有背景噪声的信号分析………………………20 六．心得体会…………………………….22 七．参考文献…………………………….23

一．摘要：

这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。

MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。

在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。

二．课程设计目的：

综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

三．设计内容：

内容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数 法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。

四．设计原理：

4.1．语音信号的采集

熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有：

a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据；

b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y；

c：wavwrite((y,fs,wavfile);创建音频文件； d：y=wavread(file);读取音频文件；

关于声音的函数还有sound()；soundsc()；等。4.2滤波器： 4．21.IIR滤波器原理

冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。

双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也 就消除了频谱混叠现象。

双线性法设计IIR数字滤波器的步骤：

1）将数字滤波器的频率指标{ k}由Wk＝（2/T）*tan(wk),转换为模拟滤波器的频率指标{k}.2)由模拟滤波器的指标设计H(s).3)由H(s)转换为H(z)21z1H(z)H(s)sT1z1

4.22.FIR滤波器原理

FIR滤波器与IIR滤波器特点不同，IIR滤波器的相位是非线性的，若需线性相位则要采用全通网络进行相位校正。而有限长单位冲激响应（FIR）数字滤波器就可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。

由于FIR系统的冲激响应就是其系统函数各次项的系数，所以设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲激响应作为H（z）的系数，冲激响应长度N就是系统函数H（z）的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。这种时域设计、频域检验的方法一般要反复几个回合，不像IIR DF设计靠解析公式一次计算成功。给出的理想滤波器频率响应是，它是w的周期函数，周期

由傅立叶反变换导出，即

hd(n)1Hd(ejw)ejwndw2，再将hd(n)与窗函数，因此可展开成傅氏级数w(n)相乘就可以得到h(n)。、的计算可采用傅氏变换的现成公式和程序，窗函数也是现成的。但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算，整个设计可用计算机编程来做。

窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有：

矩形窗

w(n)=RN(n)；

Hanning窗

；

Hamming窗

；

Blackmen窗

；

Kaiser窗。

式中Io(x)为零阶贝塞尔函数。

五．设计步骤：

5.1录制女音：

利用MATLAB中的函数录制声音。function nvyin()fs=11025;

%采样频率

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

开始录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double');

%录制声音3秒

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs);

%播放录音

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音');

%声音的存储

5.2采样语音信号并画出时域波形和频谱图

读取语音信号，画出其时域波形和频谱图，与截取后的语音信号的时域波形和频谱图比较，观察其变化。程序如下：

[x,fs,bits]=wavread('女音.wav');

%读取声音

N=length(x);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半 figure(1);subplot(2,2,1);

%把画图区域划分为2行2列，指定第一个图 plot(t, x);

%画出声音采样后的时域波形 title('原女音信号的时域波形');

%给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid;

%添加网格

y=fft(x);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,2,3);

%把画图区域划分为2行2列，指定第三个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

%添加网格

subplot(2,2,4);

%把画图区域划分为2行2列，指定第四个图 if y~=0

%判断指数是否为0

plot(k,20*log10(abs(y(n))));

%画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid;

%添加网格

%实际发出声音落后录制动作半拍的现象的解决 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]);

%截取语音信号 N=length(x1);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

figure(2);subplot(2,2,1);

%把画图区域划分为2行2列，指定第一个图 plot(t,x1);

%画出声音采样后的时域波形 title('截取后女音信号的时域波形');

%给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid;

%添加网格

y1=fft(x1);

%对信号做N点FFT变换

subplot(2,2,3);

%把画图区域划分为2行2列，指定第三个图 k=(n-1)*f0;

%频域采样点

plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

%添加网格

subplot(2,2,4);

%把画图区域划分为1行2列，指定第二个图 if y1~=0

%判断指数是否为0

plot(k,20*log10(abs(y1(n))));

%画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid;

%添加网格

原女音信号的时域波形10.5A/幅0振-0.5-10123时间/tFFT变换后声音的频谱特性FFT变换后声音的频谱特性30050A200贝/值分/幅0幅100振00200040006000-***频率/HzHz 截取后女音信号的时域波形10.5振幅/A0-0.5-10123FFT变换后声音的频谱特性50时间/tFFT变换后声音的频谱特性300200振幅/分贝幅值/A01000020004000频率/Hz6000-5002000Hz40006000

结果分析：

由原女音信号的时域波形可知录取开始时实际发出声音大概落后3500个采样点，我们把前3500点去除即可解决实际发出声音落后录制动作半拍的现象。由原女音的的频谱图和截取后声音的频谱图可看出，对声音的截取并不会影响它们频谱分布。

5.3采用双线性变换法设计IIR滤波器：

人的声音频率一般在（1~~4）kHZ之间，则我们只需要设计一个带通滤波器即可滤去声音频带以外的无用噪声，得到比较清晰的声音。根据声音的频谱图分析，设计一个带通滤波器性能指标如下：

fp1＝1000 Hz，fp2＝3000 Hz，fsc1＝500 Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000 程序如下：

%iir带通的代码： %w=2*pi*f/fs Ap=1;

%通带波纹系数

Az=100;

%最小阻带衰减

wp=[0.2 0.6];

%归一化通带数字截止频率 wz=[0.1 0.7];

%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估计契比雪夫I型滤波器阶数 [b,a]=cheby1(N,Ap,wn);

%N指定滤波器阶数，wn归一化

截 %止频率，Ap通带波动

[h,w]=freqz(b,a);

%求数字滤波器的复频率响应 figure(1);subplot(2,1,1);plot(w/pi,abs(h));

%绘制数字滤波器的频谱图 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型带通滤波器的幅频响应');subplot(2,1,2);if abs(h)~=0

%判断指数是否为0

plot(w/pi,20*log10(abs(h)));

%绘制数字滤波器频谱的分贝图 end grid;xlabel('omega/pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型带通滤波器的幅频响应');契比雪夫Ⅰ型带通滤波器的幅频响应1振幅（幅值）0.5000.10.20.50.60.70.8/契比雪夫Ⅰ型带通滤波器的幅频响应0.30.40.910振幅(分贝)-200-400-60000.10.20.30.40.5/0.60.70.80.91

5.4窗函数法设计FFR滤波器

线性相位FIR滤波器通常采用窗函数法设计。窗函数法设计FIR滤波器的基本思想是：根据给定的滤波器技术指标，选择滤波器长度N和窗函数ω（n），使其具有最窄宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h(n)。工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯泽（Kaiser）窗。

这次设计我采用的是布莱克曼来设计给定数字带通滤波器的参数如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下：

Ap=1;

%通带波纹系数 Az=100;

%最小阻带衰减 fs=10000;

%采样频率 wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));

%---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW);

%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);

%---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想带通滤波器 wdWindow=blackman(N);

%布拉克曼窗 hr=wdWindow.*hdWindow';

%点乘

n=0:N-1;

%阶数 subplot(2,2,1);stem(n,wdWindow);

%绘制布拉克曼窗时域波形 xlabel('时间');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1);

%求滤波器频率响应 subplot(2,2,3);plot(W/pi,abs(H))

%绘制滤波器频域波形 xlabel('omega/pi');ylabel('振幅');title('FIR带通滤波器幅频特性');subplot(2,2,4);

if abs(H)~=0

%判断指数是否为0

plot(W/pi,20*log10(abs(H)));

%画滤波器频谱的分贝图 end xlabel('omega/');ylabel('振幅/分贝');title('FIR带通滤波器幅频特性');grid;

%添加网格 %---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录ideallp.m function hd = ideallp(wc,N);% 理想低通滤波器的脉冲响应子程序 % hd = 点0 到 N-1之间的理想脉冲响应 % wc = 截止频率(弧度)% N = 理想滤波器的长度

tao =(N-1)/2;

% 理想脉冲响应的对称中心位置 n = [0:(N-1)];

% 设定脉冲响应长度 m = n-tao + eps;

% 加一个小数以避免零作除数

hd = sin(wc*m)./(pi*m);

% 理想脉冲响应

布拉克曼窗1振幅0.500406080时间FIR带通滤波器幅频特性500振幅/分贝20FIR带通滤波器幅频特性1.51振幅-50-100-15000.5/10.5000.5/1

5.5用IIR滤波器对信号进行滤波

用自己设计的IIR滤波器分别对采集的信号进行滤波，在Matlab中，IIR滤波器利用函数filter对信号进行滤波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半 y=fft(x);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,1);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

%iir带通的代码：

Ap=1;

%通带波纹系数

Az=100;

%最小阻带衰减

wp=[0.2 0.6];

%归一化通带数字截止频率 wz=[0.1 0.7];

%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估计契比雪夫I型滤波器阶数

[b,a]=cheby1(N,Ap,wn);

%N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x);

%对声音滤波 wavplay(x1)wavwrite(x1,'IIR滤波后女音.wav');N=length(x1);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

y=fft(x1);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,2);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

滤波前女音的频谱特性300幅值/A***030004000频率/Hz滤波后女音的频谱特性500060006040幅值/A***0频率/Hz400050006000

结果分析：

由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3400Hz的频谱成分。回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。

5.6用FIR滤波器对信号进行滤波

用自己设计的FIR滤波器分别对采集的信号进行滤波，在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波 程序如下：

[x,fs,bits]=wavread('女音.wav');N=length(x);

%计数读取信号的点数

t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

y=fft(x);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,1);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性');

%给图形加注标签说明 xlabel('频率/omega');ylabel('幅值/A');grid;

%FIR带通滤波器代码 fs=10000;wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));

%---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW);

%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);

%---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);wdWindow=blackman(N);hr=wdWindow.*hdWindow';x1=fftfilt(hr,x);

%对声音滤波 wavplay(x1)wavwrite(x1,'FIR滤波后女音.wav');N=length(x1);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

y=fft(x1);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,2);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

滤波前女音的频谱特性300200幅值/A***004000频率/l滤波后女音的频谱特性500060006040幅值/A20005001000***03000频率/Hz***0

结果分析：

由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3500Hz的频谱成分。和用IIR滤波器滤波一样，回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。5.7男女声语音信号频谱特点分析

换一个男音录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点。程序如下：

[x,fs,bits]=wavread('女音.wav');

%读取声音

N=length(x);

%计数读取信号的点数 t=(1:N)/fs;

f0=fs/N;

n=1:N/2;

y=fft(x);

k=(n-1)*f0;

subplot(2,1,1);

plot(k,abs(y(n)));

title('FFT变换后女音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid;

[x,fs,bits]=wavread('明明.wav');

N=length(x);

t=(1:N)/fs;

f0=fs/N;

n=1:N/2;

y=fft(x);

k=(n-1)*f0;

subplot(2,1,2);

plot(k,abs(y(n)));

title('FFT变换后男音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid;

%信号的时域采样点

%采样间隔

%取信号的一半

%对信号做N点FFT变换

%频域采样点

%把画图区域划分为2行1列，指定第一个图%绘制原始语音信号的幅频响应图

%给图形加注标签说明

%添加网格

%读取声音

%计数读取信号的点数

%信号的时域采样点

%采样间隔

%取信号的一半

%对信号做N点FFT变换

%频域采样点

%把画图区域划分为2行1列，指定第二个图%绘制原始语音信号的幅频响应图

%给图形加注标签说明

%添加网格

axis([0 6000 0 300])；

%改变横纵坐标便于比较频谱图

FFT变换后女音的频谱特性300200幅值/A***00频率/FFT变换后男音的频谱特性***200幅值/A***00频率/400050006000

结果分析：

通过比较上面女音频谱图和男音频谱图可知，男音的频谱集中在低频部分，高频成分底，谱线较平滑，声音听起来低沉。5.8有背景噪声的信号分析

从硬盘中把一段噪声（频谱能量集中在某个小范围内）叠加到语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除； 程序如下：

z=wavread('女音.wav',[1 24000]);

%读取声音在1-24000之间 f=wavread('noise.wav',[1 24000]);x=z+f;wavplay(x);fs=11025;N=length(x);f0=fs/N;

%采样间隔

n=1:N;

%取信号的一半 y=fft(x,N);%对信号做N点FFT变换

k=(n-1)*f0;

%频域采样点

subplot(2,1,1);

%把画图区域划分为1行2列，指定第二个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('加噪声后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格

%iir带通滤波器的代码：

Ap=1;

%通带波纹系数

Az=70;

%最小阻带衰减

wp=[0.2 0.7];

%归一化通带数字截止频率 wz=[0.1 0.8];

%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估计契比雪夫I型滤波器阶数

[b,a]=cheby1(N,Ap,wn);

%N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x);

%对声音滤波 wavplay(x1)；

wavwrite(x1,'滤除噪音后女音.wav');N=length(x1);f0=fs/N;

%采样间隔 n=1:N;

%取信号的一半

y1=fft(x1,N);

%对信号做fs点FFT变换

subplot(2,1,2);

%把画图区域划分为1行2列，指定第二个图 k=(n-1)*f0;

%频域采样点

plot(k,abs(y1(n)));

%绘制原始语音信号的幅频响应图 title('滤除噪声后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格

加噪声后声音的频谱特性3000幅值/A***0008000频率/Hz滤除噪声后声音的频谱特性***030幅值/A***000频率/Hz80001000012000

结果分析

观察加噪声后声音的频谱图可知，噪音频率主要在4000Hz处，只要我们设计一个，滤波器滤除大概在4000Hz的频谱即可，回放滤波后的语音信号，可证噪音基本滤除。

六．心得体会：

通过这次课程设计，让我对MATLAB的基本应用有了更深的了解，还有数字信号处理在MATLAB中的一些函数的用法。通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

在这次实验中，也遇到了很多问题，比如画信号频谱的分贝图时（20*log10(abs(y))）指数为零时的处理。滤波器的设计也花了好大的功夫，刚开始不会设计参数，一头雾水，通过同学的指导和讨论，得知通过观察信号的频谱图，看噪音频率集中在那一部分，设计滤波器把其滤除即可。可反复设置参数直到滤波后语音信号的效果好为止。

七．参考文献：

（1）《MATLAB LabVIEW SystemView》翁剑枫 叶志前 编著, 机械工业出版社；

（2）《MATLAB及在电子信息课程中的应用》陈怀琛 吴大正 高西全编著，电子工业出版社；

（3）《MATLAB在数字信号处理中的应用》（弟2版）薛年喜 编著，清华大学出版社；

（4）《MATLAB扩展编程》何强 何英

编著，清华大学出版社；（5）《MATLAB7简明教程》吴清 曹辉林 编著，清华大学出版社；（6）MATLAB5.3精要.编程及高级应用》程卫国 冯峰 王雪梅 刘艺 编著，机械工程出版社。

第二篇：数字信号处理课后习题Matlab作业

数字信号处理MATLAB

第1页

习题数字信号处理MATLAB习题

M1-1 已知g1(t)cos(6t)，g2(t)cos(14t)，g3(t)cos(26t)，以抽样频率fsam10Hz对上述三个信号进行抽样。在同一张图上画出g1(t)，g2(t)和g3(t)及抽样点，对所得结果进行讨论。

解：

第2页

从以上两幅图中均可看出，三个余弦函数的周期虽然不同，但它们抽样后相应抽样点所对应的值都相同。那么这样还原回原先的函数就变成相同的，实际上是不一样的。这是抽样频率太小的原因，我们应该增大抽样频率才能真实还原。如下图：f=50Hz

第3页

程序代码

f=10;

t=-0.2:0.001:0.2;g1=cos(6.*pi.*t);g2=cos(14.*pi.*t);g3=cos(26.*pi.*t);k=-0.2:1/f:0.2;h1=cos(6.*pi.*k);h2=cos(14.*pi.*k);h3=cos(26.*pi.*k);% subplot(3,1,1);

% plot(k,h1,'r.',t,g1,'r');% xlabel('t');% ylabel('g1(t)');% subplot(3,1,2);

% plot(k,h2,'g.',t,g2,'g');% xlabel('t');% ylabel('g2(t)');% subplot(3,1,3);

% plot(k,h3,'b.',t,g3,'b');% xlabel('t');% ylabel('g3(t)');

plot(t,g1,'r',t,g2,'g',t,g3,'b',k,h1,'r.',k,h2,'g.',k,h3,'b.')

第4页

xlabel('t');ylabel('g(t)');

legend('g1(t)','g2(t)','g3(t)');

M2-1 利用DFT的性质，编写一MATLAB程序，计算下列序列的循环卷积。

(1)g[k]={1,-3,4,2,0,-2,},h[k]={3,0,1,-1,2,1}；(2)x[k]=cos(k/2),y[k]=3k,k=0,1,2,3,4,5。解：(1)循环卷积结果

6.0000-3.0000 17.0000-2.0000 7.0000-13.0000

程序代码

第5页

g=[1-3 4 2 0-2];h=[3 0 1-1 2 1];l=length(g);L=2*l-1;GE=fft(g,L);HE=fft(h,L);y1=ifft(GE.*HE);for n=1:l

if n+l<=L

y2(n)=y1(n)+y1(n+l);else

y2(n)=y1(n);

end end y2

stem(0:l-1,y2)xlabel('k')ylabel('y(k)')title('循环卷积')

(2)循环卷积结果

-71.0000-213.0000 89.0000 267.0000 73.0000 219.0000

第6页

程序代码

k=0:5;

x=cos(pi.*k./2);y=3.^k;l=length(x);L=2*l-1;GE=fft(x,L);HE=fft(y,L);y1=ifft(GE.*HE);for n=1:l

if n+l<=L

y2(n)=y1(n)+y1(n+l);

else

y2(n)=y1(n);

end end y2

stem(0:l-1,y2)xlabel('k')ylabel('y’(k)')title('循环卷积')

第7页

M2-2 已知序列x[k]cos(k/2N),|k|N

0,其他(1)计算序列DTFT的表达式X(ej),并画出N=10时，X(ej)的曲线。

(2)编写一MATLAB程序，利用fft函数，计算N=10时，序列x[k]的DTFT在m2m/N的抽样值。利用hold函数，将抽样点画在X(ej)的曲线上。

解：

(1)X(e)DTFT{x[k]}jkx[k]ejkkNcos(k/2N)eNjk

程序代码

N=10;k=-N:N;

x=cos(k.*pi./(2*N));W=linspace(-pi,pi,512);

第8页

X=zeros(1,length(W));for k=-N:N

X1=x(k+N+1).*exp(-j.*W.*k);X=X+X1;end

plot(W,abs(X))xlabel('W');ylabel('abs(X)');

(2)

程序代码

N=10;k=-N:N;

x=cos(k.*pi./(2*N));X_21=fft(x,21);L=-10:10;

W=linspace(-pi,pi,1024);X=zeros(1,length(W));for k=-N:N

X1=x(k+N+1).*exp(-j.*W.*k);X=X+X1;end

第9页

plot(W,abs(X));hold on;

plot(2*pi*L/21,fftshift(abs(X_21)),'o');xlabel('W');ylabel('abs(X)');

M2-3 已知一离散序列为x[k]Acos0kBcos[(0)k]。用长度N=64的Hamming窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用Hamming窗能分辨的最小的谱峰间隔wc的值。

解：f1=100Hz f2=120Hz时

2中cN

f2=140Hz时

第10页

f2=160Hz时

第11页

由以上三幅图可见

f2=140Hz时，各谱峰可分辨。则f又

wc2N

40Hz

且

wT2fT2401 800所以c=3.2（近似值）

程序代码

N=64;L=1024;

f1=100;f2=160;;fs=800;

A=1;B1=1;B2=0.5;B3=0.25;B4=0.05;T=1/fs;ws=2*pi*fs;k=0:N-1;

x1=A*cos(2*pi*f1*T*k)+B1*cos(2*pi*f2*T*k);x2=A*cos(2*pi*f1*T*k)+B2*cos(2*pi*f2*T*k);x3=A*cos(2*pi*f1*T*k)+B3*cos(2*pi*f2*T*k);x4=A*cos(2*pi*f1*T*k)+B4*cos(2*pi*f2*T*k);hf=(hamming(N))';x1=x1.*hf;x2=x2.*hf;x3=x3.*hf;x4=x4.*hf;

X1=fftshift(fft(x1,L));X2=fftshift(fft(x2,L));X3=fftshift(fft(x3,L));X4=fftshift(fft(x4,L));

W=T*(-ws/2+(0:L-1)*ws/L)/(2*pi);subplot(2,2,1);plot(W,abs(X1));title('A=1,B=1');xlabel('W');ylabel('X1');subplot(2,2,2);

第12页

plot(W,abs(X2));title('A=1,B=0.5');xlabel('W');ylabel('X2');subplot(2,2,3);plot(W,abs(X3));title('A=1,B=0.25');xlabel('W');ylabel('X3');subplot(2,2,4);plot(W,abs(X4));title('A=1,B=0.05');xlabel('W');ylabel('X4');

M2-4 已知一离散序列为x[k]cos0k0.75cos1k,0k63。其中, 02/15，12.3/15。

(1)对x[k]做64点FFT, 画出此时信号的谱。

(2)如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补0而分辨出两个谱峰。通过编程进行证实，并解释其原因。

解：(1)

第13页

程序代码

W0=2*pi/15;W1=2.3*pi/15;N=64;k=0:N-1;

x=cos(W0*k)+0.75*cos(W1*k);X=fft(x);

plot(k/N,abs(X));grid on;

title('64点FFT');

(2)

第14页

第15页

由以上三幅图看出：不能对(1)中的64点信号补零而分辨出两个谱峰，这样的方法只能改变屏幕分辨率，但可以通过加hamming窗来实现对谱峰的分辨。程序代码

W0=2*pi/15;W1=2.3*pi/15;N=64;L=1024;k=0:N-1;

x=cos(W0*k)+0.75*cos(W1*k);X=fft(x,L);

plot((0:L-1)/N,abs(X));grid on;

title('1024点FFT');

M2-5 已知一连续信号为x(t)=exp(-3t)u(t)，试利用DFT近似分析

第16页

其频谱。若要求频率分辨率为1Hz，试确定抽样频率fsam、抽样点数N以及持续时间Tp。

解：

本题使用矩形窗，则Nfsamfsam1fsam，Tp1 f1f

第17页

由以上三幅图可以看出当fsam越来越大时，近似值越来越接近

第18页

于实际值。即fsam越大拟合效果越好，造成的混叠也是在可以允许的范围内。程序代码

fs=100;ws=2*pi*fs;Ts=1/fs;N=fs;

x=exp(-3*Ts*(0:N-1));y=fft(x,N);l=length(y);

k=linspace(-ws/2,ws/2,l);

plot(k,Ts*fftshift(abs(y)),'b:');hold on;

w=linspace(-ws/2,ws/2,1024);y1=sqrt(1./(9+w.^2));plot(w,y1,'r')

title('fs=100Hz时的频谱')legend('近似值','实际值);

M2-6 试用DFT近似计算高斯信号g(t)exp(dt2)的频谱抽样值。

π2通过和频谱的理论值G(j)exp()比较，讨论如何根据时域的信

d4d号来恰当地选取截短长度和抽样频率使计算误差能满足精度要求。

解：

第19页

第20页

由以上三幅图可以看出：

当时域截取长度相同时，抽样间隔越小时误差越小，当抽样间隔一定时，时域截取长度越长，误差越小。当取抽样间隔为1S，时域截取长度为2S时，误差较大，绝对误差在0.5左右；当抽样间隔为0,5S，时域截取长度为2S时，误差比间隔为1S时小，绝对误差不大于0.2；当抽样间隔为0.5S时域截取长度为4S时，误差更小，绝对误差不大于0.04。因为时域截取长度越长，保留下来的原信号中的信息越多，抽样间隔越小，频谱越不容易发生混叠，所以所得频谱与理论值相比，误差更小。

程序代码

Ts=0.5;N=4;N0=64;

k=(-N/2:(N/2))*Ts;

第21页

x=exp(-pi*(k).^2);X=Ts*fftshift(fft(x,N0));

w=-pi/Ts:2*pi/N0/Ts:(pi-2*pi/N0)/Ts;XT=(pi/pi)^0.5*exp(-w.^2/4/pi);subplot(2,1,1)

plot(w/pi,abs(X),'-o',w/pi,XT);xlabel('omega/pi');ylabel('X(jomega)');

legend('试验值','理论值');

title(['Ts=',num2str(Ts)subplot(2,1,2)plot(w/pi,abs(X)-XT)ylabel('实验误差')

xlabel('omega/pi');

'N=',num2str(N)]);第22页

' '

第三篇：Matlab的数字信号处理课程实验设计的论文

摘要:本文设计了一个基于Matlab的“数字信号处理”课程综合性实验。该实验把“数字信号处理”课程中的许多离散的知识点串接了起来，包括采样、量化、滤波器设计、滤波器实现、DFT/FFT和滤波器的有限字长效应等。教学实践表明该实验有利于巩固学生课堂上学到的理论知识，提高学生的理论联系实际的能力和动手解决问题的能力。

关键词:数字信号处理;综合性实验;Matlab

0引言

“数字信号处理”课程的主要内容包括z变换、离散傅里叶变换(DFT)、快速傅里叶变换(FFT)、数字滤波器设计和实现以及数字信号处理中的有限字长效应等等［1］。在学习理论知识的同时或之后，引入实验将有助于学生更好地理解和掌握课程内容［2－3］。笔者在教学过程中，设计了Matlab综合性实验。该实验在不失趣味性的同时，能把该课程中许多分散的知识点串接起来。教学实践表明，该实验可以帮助学生更深入地理解本门课程，取得了较好的教学效果。

1综合实验内容设计

笔者所设计的Matlab实验如下:对下式所示的输入信号进行滤波。x=sin(100πt)+sin(480πt)(1)具体步骤为(1)将输入的模拟信号x进行采样和量化，得到12位精度的数字信号;(2)设计一个低通无限冲激响应(IIＲ)滤波器，将输入信号中的240Hz的干扰滤除，要求滤波器的输出信号中240Hz处的噪声功率比50Hz处的信号功率低60dB。(3)设计一个高通有限冲激响应(FIＲ)滤波器，将输入信号中的50Hz的干扰滤除，要求滤波器的输出信号中50Hz处的噪声功率比240Hz处的信号功率低60dB。(4)对于上述两个滤波器，要求:给出理想滤波器的传输函数及频率响应;给出系数量化后所得的新的滤波器的传输函数及频率响应;确定滤波器实现所采用的结构，并给出该结构中所用加法器和乘法器的位数;将输入的数字信号通过前一步实现的滤波器，画出输出信号的频谱，确保滤波器性能满足设计要求。顺利完成上述Matlab实验，需要解决以下问题:(1)采样频率和FFT点数的选取:根据采样定理，采样频率只要不低于信号中所包含的最高频率的两倍，就可以从采样后的离散时间信号中恢复出原始的模拟信号。根据式(1)，采样频率只要不小于480Hz即可。但是当需要使用FFT对信号进行频谱分析时，在确定采样频率时，除了要满足采样定理外，还需要考虑其他条件。例如:在做FFT时，信号频率应为频率分辨率的整数倍，这样才能准确地从频谱中看到该频率信号的功率，避免谱泄漏，即下式中的k应为整数:k=ffs=N(2)其中f，fs和N分别为信号频率、采样频率和FFT的点数。fs/N为频率分辨率，N一般为2的幂次方。在k不为整数时，为了减小谱泄漏的影响，可以在做FFT之前对采样所得的信号进行加窗处理［1］。(2)模数转换器的实现:实验中要求对输入信号进行量化，得到12位精度的数字信号。在将输入信号进行量化时，涉及到如何确定模数转换器的满量程范围、结构、量化方式(舍入还是截断)以及如何进行有符号数的量化等。(3)IIＲ滤波器类型的选择和设计:双线性变换是设计数字IIＲ滤波器的常用方法。它首先要将所要设计的数字滤波器的归一化边界角频率进行预畸变，然后再设计出满足性能要求的模拟滤波器。模拟滤波器有四种类型，分别为巴特沃斯滤波器，切比雪夫I型滤波器、切比雪夫II型滤波器以及椭圆滤波器。只有了解了这四种滤波器的特性，才能根据实际需求来选择合适的滤波器类型。在选择好滤波器类型后，将滤波器的性能指标输入相应的Matlab函数，就可以得到滤波器的传输函数，完成滤波器的设计。以椭圆滤波器为例，可以依次调用函数elli-pord()，函数ellipap()和函数zp2tf()来获得滤波器的阶数、零极点、增益和s域传输函数;也可以直接调用函数ellip()来得到滤波器的s域传输函数。最后再通过调用函数bilinear()得到相应数字滤波器的传输函数。(4)FIＲ滤波器的设计:在用窗函数法来设计FIＲ滤波器时，首先要根据滤波器的性能参数(如过渡带宽度、阻带衰减等)选取合适的窗函数以及确定窗函数的长度，之后将得到的窗函数与理想滤波器的单位脉冲响应序列相乘得到FIＲ滤波器的单位脉冲响应序列。以Kaiser窗为例，在Matlab中，函数kaiserord()用于预估FIＲ滤波器的阶数，函数kaiser()用于产生相应长度的Kaiser窗函数，函数fir1()用于实现采用该Kaiser窗设计的FIＲ滤波器，输出为滤波器的单位脉冲响应序列。(5)滤波器的实现:在用硬件实现滤波器时，必须考虑滤波器的有限字长效应，即滤波器系数的量化、滤波器中加法器和乘法器的有限字长效应以及运算结果的有限字长等等。滤波器的实现结构有直接型、级联型和并联型等。由于IIＲ滤波器存在量化噪声的积累，所以在选择结构时，需要考虑各种结构对有限字长效应的灵敏度。高阶IIＲ滤波器通常采用级联型或并联型结构来实现。Matlab中的函数residuez(B，A)用于计算传输函数B(z)/A(z)的留数、极点和直接项，从而得到有理式的部分分式展开;利用传输函数的部分分式展开，并通过适当的合并，可以得到滤波器的并联型结构。函数tf2sos()则可用于将传输函数转换成二阶节，得到滤波器的级联型结构。图3给出了系数量化前后高通滤波器的频率响应。为了能够判断所设计和实现的滤波器的性能是否达到设计指标，需要对滤波器的输出序列做N点的FFT。这时需要注意两点:一要能正确地区分输出序列中的暂态响应部分和稳态响应部分;二要从稳态响应部分选取连续的N个输出值做N点的FFT。

2教学反馈

根据学生上交的实验报告，从他们所写的实验收获和实验心得可以看出这个实验对他们学好这门功课所起的作用。总结如下:(1)本次实验是FIＲ滤波器与IIＲ滤波器的设计，综合使用了大量数字滤波器的设计方法，比如双线性变换法，窗函数法等，加深了对课堂学习的理论知识的理解，如IIＲ和FIＲ滤波器的优缺点、滤波器的暂态响应和稳态响应、各种模拟滤波器的性能比较以及各种窗函数之间的差异等。(2)学生对采样定理和FFT有了更深的认识，明白了采样频率、FFT点数等对频谱分析结果的影响，并通过不断的摸索与尝试，总结出了使用FFT时的一些注意事项。(3)对数字信号处理中的有限字长效应有了更加直观的体会，认识到在设计滤波器的传输函数时，需要考虑量化对滤波器性能的影响，设计指标需要留出一定的裕量。(4)提高了用Matlab实现数字信号处理功能的能力，包括:熟悉了使用Matlab设计FIＲ和IIＲ滤波器的流程;学会使用Matlab中的一些函数，如fft，cheb1ord，cheby，bilinear，fir1等;学会了用Matlab编写程序来实现指定结构的滤波器;学会了从时域和频域观察滤波器的输出是否正确以及是否达到性能要求等。总而言之，通过这次实验，使学生真正了解了如何利用Matlab来进行滤波器的设计，感觉受益匪浅，对他们学好“数字信号处理”课程很有帮助。

3结语

笔者所设计的基于Matlab的综合性实验涵盖了“数字信号处理”课程中的主要知识点。从学生反馈的意见可以看出，本实验取得了良好的教学效果，这有利于提高学生学习兴趣以及增强他们解决实际问题的能力。

参考文献:

［1］程佩青，数字信号处理教程［M］，北京:清华大学出版社，2007．

［2］曹建玲，刘焕淋，雷宏江．基于MATLAB的“数字信号处理”仿真实验［J］．北京:中国电力教育，2012(32):88－89．

［3］易婷．“数字信号处理”课程课内配套实验的设计［J］．南京:电气电子教学学报，2013，35(4):89－90．

第四篇：《数字信号处理原理及实现》课程小结

时间过得好快，转眼半学期结束了。这半学期数字信号的学习让我受益匪浅。前两章和信号与线性系统相关，介绍了离散时间信号与系统的时域分析方法最深刻的是采样，时域采样定理与采样恢复和采样内插公式。第二章介绍了离散时间信号与系统的频域分析，DTFT与Z变换，系统函数的零极点分布。第三章主要讲了离散傅里叶变换DFT及其性质，和频域采样定理。第四章介绍了傅里叶快速变换FFT，熟悉了其原理特点及方法。五六两章分别介绍了IIR和FIR滤波器，知道了IIR的脉冲响应不变法与双线性变换法及其优缺点，并学会其MATLAB应用设计滤波器，FIR的窗函数法与频域采样法设计滤波器及其MATLAB实现。第七章主要介绍了IIR和FIR滤波器的基本网络结构，通过老师上课习题的练习基本掌握了其结构图的画法。

先说说对课程的建议吧，张晓光老师是个很负责讲课思路也很清晰的老师，知道从学生的角度来讲问题，根据学生的反应来调整课程进度。我们都很喜欢这样的老师，老师开新课之前总是列提纲复习上节课讲的知识，每章结束都根据章节的重要性开一节总结课，这种方式个人觉得很好，希望老师坚持。但是，感觉老师讲题讲的不是很多，或许是课时原因，但我觉得每章结束后开一节例题课，把知识点融进去，毕竟大学生现在做题比较少，这样强制一下效果会更好。这次考试的试题觉得有不少都见过，有的是课后题，但做起来还是有点吃力，应该就是习题练的少，计算跟不上去。至于教材，我觉得编的很好，每章都有相关的MATLAB编程方法，在原理讲清之后就来实践，免去了学生盲目做实验，提高了效率。还有就是老师也很重视实验，总是把相关的MATLAB语句语义讲解清楚，这样我们在编程序时也就相对容易点。但好像老师讲程序时都注重程序的意思了，希望老师以后再讲程序时把它先部分后整体，就是在讲完程序意思后把程序设计思路或框架结构，及各部分要实现什么再讲讲，这样有助于学生设计时设计思路更清晰。再说说考试，老师分卷面成绩和实验成绩及平时成绩，将实验单独考试，可见对实验的重视，也说明MATLAB的重要性，这样确实提高了学生的重视心理，虽然实验做完了，但做完50道题并看完相关讲解，我又收获了不少，理清了设计方法与思路，所以我觉的考试方式还是挺不错的，锻炼了我们各方面的知识。

数字信号课程结束了，真希望您还能教我们别的课。

小组成员：陈文斌、李亚伟、王猛、汪子雄、吴官宝

第五篇：数字信号处理实验-FFT的实现

实

验

报

告

学生姓名：

学 号：

指导教师：

一、实验室名称：数字信号处理实验室

二、实验项目名称：FFT的实现

三、实验原理：

一．FFT算法思想：

1．DFT的定义：

对于有限长离散数字信号{x[n]}，0  n  N-1，其离散谱{x[k]}可以由离散付氏变换（DFT）求得。DFT的定义为：

N1X[k]通常令ej2Nx[n]en0j2Nnk，k=0,1,…N-1 WN，称为旋转因子。

2．直接计算DFT的问题及FFT的基本思想：

由DFT的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（N-1）次加法。因此，对于一些相当大的N值（如1024）来说，直接计算它的DFT所作的计算量是很大的。

FFT的基本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。例如，若N为偶数，将原有的N

22点序列分成两个（N/2）点序列，那么计算N点DFT将只需要约[(N/2)·2]=N/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。

3．基2按时间抽取（DIT）的FFT算法思想：

设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

将长度为N2L的序列x[n](n0,1,...,N1)，先按n的奇偶分成两组：

x[2r]x1[r]x[2r1]x2[r]，r=0,1,…,N/2-1 DFT化为：

N1N/21N/21X[k]DFT{x[n]}N/21n0x[n]WnkN2rkr0x[2r]W2rkNr0x[2r1]WN(2r1)kN/21r0N/21x1[r]Wx1[r]W2rkNWWkNr0N/21x2[r]WN

r0rkN/2kNr0x2[r]WN/22rkrk上式中利用了旋转因子的可约性，即：WNN/21NrkN/21rkWN/2。又令

rkX1[k]r0x[1r]W,/X2[k]2r0x[r]WN2，则上式可以写成： /2X[k]X1[k]WNX2[k](k=0,1,…,N/2-1)

k可以看出，X1[k],X2[k]分别为从X[k]中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了X[k]前N/2点的DFT值，还需要继续利用X1[k],X2[k]表示X[k]的后半段本算法推导才完整。利用旋转因子的周期性，有：WN/2WN/2X1[N2N/21rkr(kN/2)，则后半段的DFT值表达式：

rkk]r0x1[r]W2N/2r(Nk)N/21r0x1[r]WN/2X1[k]，同样，X2[N2k]X2[k]

（k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表达式获得，中间还利用到WN(N2k)NWN2Wk得到后半段的X[k]值表达式W，k为：X[k]X1[k]WNkX2[k]（k=0,1,…,N/2-1）。

这样，通过计算两个N/2点序列x1[n],x2[n]的N/2点DFTX1[k],X2[k]，可以组合得到N点序列的DFT值X[k]，其组合过程如下图所示：

X1[k] X1[k]WNkX2[k]

X2[k] WNnk-1 X1[k]WNkX2[k]

比如，一个N = 8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：

x(0)W0x(1)W0x(2)W0x(3)W2W0W1W0x(5)W0x(6)W0x(7)W2X(7)W3X(6)W2X(5)X(3)X(2)X(1)X(0)x(4)X(4)

4．基2按频率抽取（DIF）的FFT算法思想：

设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

在把X[k]按k的奇偶分组之前，把输入按n的顺序分成前后两半：

N1N/21nkNN1X[k]DFT{x[n]}N/21N/21x[n]Wn0(nn0N2)kx[n]WnkNnN/2x[n]WNnkn0N/21x[n]WnkNn0x[nNkN2]WNnk

Nn0[x[n]x[nN2NkN2]W2N]WN,k0,1,...,N1因为W2N1，则有WX[k](1)，所以：

kkN/21n0[x[n](1)x[nN2]]WN,k0,1,...,N1

nk按k的奇偶来讨论，k为偶数时：

N/21X[2r]n0[x[n]x[nN2]]WN,k0,1,...,N1 N22rnN/21k为奇数时：X[2r1]前面已经推导过WNN/21n0[x[n]x[n]]WN(2r1)n,k0,1,...,N1

2rkWN/2，所以上面的两个等式可以写为：

N2]]WN/2,r0,1,...,N/21 N2rnrkX[2r]n0[x[n]x[nN/21X[2r1]n0{[x[n]x[n]]WN}WN/2,r0,1,...,N/21

nnr通过上面的推导，X[k]的偶数点值X[2r]和奇数点值X[2r1]分别可以由组合而成的N/2点的序列来求得，其中偶数点值X[2r]为输入x[n]的前半段和后半段之和序列的N/2点DFT值，奇数点值X[2r1]为输入x[n]的前半段和后半段之差再与WN相乘序列的N/2点DFT值。

令x1[n]x[n]x[nN/21nN2]，x2[n][x[n]x[nN/21N2]]WN，则有：

nX[2r]n0x1[n]WrnN/2,X[2r1]n0x2[n]WrnN/2,r0,1,...,N21

这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示：

x[n] x[n]x[nN2]

x[nN2]-1 WNn [x[n]x[nN2]]WNn

二．在FFT计算中使用到的MATLAB命令：

函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R

四、实验目的：

离散傅氏变换（DFT）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。

本实验通过直接计算DFT，利用FFT算法思想计算DFT，以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。

五、实验内容：

a)计算实数序列x(n)cos516n,0n256的256点DFT。

b)计算周期为1kHz的方波序列（占空比为50％，幅度取为＋/-512，采样频率为25kHz，取256点长度）256点DFT。

六、实验器材（设备、元器件）：

安装MATLAB软件的PC机一台，DSP实验演示系统一套。

七、实验步骤：

（1）先利用DFT定义式，编程直接计算2个要求序列的DFT值。

（2）利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值。（3）(拓展要求)不改变序列的点数，仅改变DFT计算点数（如变为计算1024点DFT值），观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。

（4）利用FFT的基本思想（基2－DIT或基2－DIF），自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果进行对比。

（5）(拓展要求)尝试对其他快速傅立叶变换算法（如Goertzel算法）进行MATLAB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行对比。

（6）（拓展要求）在提供的DSP实验板上演示要求的2种序列的FFT算法（基2－DIT），用示波器观察实际计算出来的频谱结果，并与理论结果对比。

八、实验数据及结果分析：

程序：（1）对要求的2种序列直接进行DFT计算的程序

（2）对要求的2种序列进行基2－DIT和基2－DIF FFT算法程序（3）对要求的2种序列用MATLAB中提供的FFT函数进行计算的程序

结果：（1）对2种要求的序列直接进行DFT计算的频域波形

（2）对2种要求的序列进行基2－DIT和基2－DIF FFT算法频域波形（3）对2种要求的序列用MATLAB中提供的FFT函数计算的频域波形。（4）（拓展要求）分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。

（5）（拓展要求）保持序列点数不变，改变DFT计算点数（变为1024点），观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。

九、实验结论：

十、总结及心得体会：

十一、对本实验过程及方法、手段的改进建议：