DSP研究性学习报告.docx
- 文档编号:6451751
- 上传时间:2023-01-06
- 格式:DOCX
- 页数:17
- 大小:83.26KB
DSP研究性学习报告.docx
《DSP研究性学习报告.docx》由会员分享,可在线阅读,更多相关《DSP研究性学习报告.docx(17页珍藏版)》请在冰豆网上搜索。
DSP研究性学习报告
《数字信号处理》课程研究性学习报告
姓名王佳
学号09211013
同组成员李禹霏
厚福佳
胡远洋
白雪尧
辛未
指导教师黄琳琳
时间2011年4月9日
DFT近似计算信号频谱专题研讨
【目的】
(1)掌握利用DFT近似计算不同类型信号频谱的原理和方法。
(2)理解误差产生的原因及减小误差的方法。
(3)培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】
1.利用DFT分析x(t)=Acos(2pf1t)+Bcos(2pf2t)的频谱,其中f1=100Hz,f2=120Hz。
(1)A=B=1;
(2)A=1,B=0.2
要求选择不同的DFT参数及窗函数,并对实验结果进行比较,总结出选择合适DFT参数的原则.
【题目分析】 频率分辨率在信号处理领域是否达到要求,会直接影响分析的结果。
在最小采样率的条件下,可以用补零DFT来提高信号可视频率分辨率,但信号真正的频率分辨率并没有得到改善。
通过对补零DFT和信号采样点数改变的分析研究,可以区分可视分辨率和真正频率分辨率两个重要的概念
【结果分析】对于窗函数的选择,应考虑被分析信号的性质与处理要求。
如果仅要求精确读出主瓣频率,而不考虑幅值精度,则可选用主瓣宽度比较窄而便于分辨的矩形窗,例如测量物体的自振频率等;如果分析窄带信号,且有较强的干扰噪声,则应选用旁瓣幅度小的窗函数,如汉宁窗、三角窗等;对于随时间按指数衰减的函数,可采用指数窗来提高信噪比。
【自主学习内容】不同的窗函数对信号频谱的影响是不一样的,这主要是因为不同的窗函数,产生泄漏的大小不一样,频率分辨能力也不一样。
信号的截断产生了能量泄漏,而用FFT算法计算频谱又产生了栅栏效应,从原理上讲这两种误差都是不能消除的,但是我们可以通过选择不同的窗函数对它们的影响进行抑制。
(矩形窗主瓣窄,旁瓣大,频率识别精度最高,幅值识别精度最低;布莱克曼窗主瓣宽,旁瓣小,频率识别精度最低,但幅值识别精度最高)
【阅读文献】
【发现问题】
【问题探究】
【仿真程序】
N=30;
L=512;
f1=100;
f2=120;
fs=600;
T=1/fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
y=fft(x,L);
mag=abs(y);
f=(0:
length(y)-1)*fs/length(y);
plot(f(1:
L/2),mag(1:
L/2));
N=600;
L=512;
f1=100;
f2=120;
fs=600;
T=1/fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
y=fft(x,L);
mag=abs(y);
f=(0:
length(y)-1)*fs/length(y);
plot(f(1:
L/2),mag(1:
L/2));
N=20;
L=512;
f1=100;
f2=120;
fs=600;
T=1/fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
y=fft(x,L);
mag=abs(y);
f=(0:
length(y)-1)*fs/length(y);
plot(f(1:
L/2),mag(1:
L/2));
A=1;B=1用哈明窗
N=30;%数据的长度
L=512;%DFT的点数
f1=100;
f2=120;
fs=600;%抽样频率
T=1/fs;%抽样间隔
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
wh=(hamming(N));
X=x*wh;
X=fftshift(fft(x,L));
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X))
使用矩形窗
N=30;
L=512;
f1=100;
f2=120;
fs=600;
T=1/fs;
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
wh=(boxcar(N));
X=x*wh;
X=fftshift(fft(x,L));
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X))
A=1,B=0.2
矩形窗
N=30;%Êý¾ÝµÄ³¤¶È
L=512;%DFTµÄµãÊý
f1=100;
f2=120;
fs=600;%³éÑùƵÂÊ
T=1/fs;%³éÑù¼ä¸ô
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);
wh=(boxcar(N));
X=x*wh;
X=fftshift(fft(x,L));
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X))
用哈明窗
N=30;%Êý¾ÝµÄ³¤¶È
L=512;%DFTµÄµãÊý
f1=100;
f2=120;
fs=600;%³éÑùƵÂÊ
T=1/fs;%³éÑù¼ä¸ô
ws=2*pi*fs;
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);
wh=(hamming(N));
X=x*wh;
X=fftshift(fft(x,L));
w=(-ws/2+(0:
L-1)*ws/L)/(2*pi);
plot(w,abs(X))
【研讨题目】
2.分析矩形窗、Hann窗、Hamming窗、Blackman窗、Kaiser窗的频谱,并进行比较。
【题目分析】
矩形窗
矩形窗属于时间变量的零次幂窗。
矩形窗使用最多,习惯上不加窗就是使信号通过了矩形窗。
这种窗的优点是主瓣比较集中,缺点是旁瓣较高,并有负旁瓣,导致变换中带进了高频干扰和泄漏,甚至出现负谱现象。
Hanning窗
汉宁窗又称升余弦窗,汉宁窗可以看作是3个矩形时间窗的频谱之和,或者说是3个sinc(t)Hamming窗型函数之和,而括号中的两项相对于第一个谱窗向左、右各移动了π/T,从而使旁瓣互相抵消,消去高频干扰和漏能。
可以看出,汉宁窗主瓣加宽并降低,旁瓣则显著减小,从减小泄漏观点出发,汉宁窗优于矩形窗.但汉宁窗主瓣加宽,相当于分析带宽加宽,频率分辨力下降
Hamming窗
海明窗也是余弦窗的一种,又称改进的升余弦窗。
海明窗与汉宁窗都是余弦窗,只是加权系数不同。
海明窗加权的系数能使旁瓣达到更凯泽(Kaiser)窗
小。
分析表明,海明窗的第一旁瓣衰减为一42dB.海明窗的频谱也是由3个矩形时窗的频谱合成,但其旁瓣衰减速度为20dB/(10oct),这比汉宁窗衰减速度慢。
海明窗与汉宁窗都是很有用的窗函数。
Kaiser窗是一种可以任意调整窗谱主瓣宽度与旁瓣衰减比例的函数。
以瞬时无功功率理论为基础,提出一种谐波电流的实时检测新方法,用Kaiser滑动时窗截取谐波电流信号,通过对窗函数参数的选定,提高了谐波电流检测的准确性。
1.矩形窗:
利用w=boxcar(n)的形式得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,它的元素由窗函数的值组成。
‘w=boxcar(n)’等价于‘w=ones(1,n)’.
2.三角窗:
利用w=triang(n)的形式得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,它的元素由窗函数的值组成。
w=triang(N-2)等价于bartlett(N)。
3.汉宁窗:
利用w=hanning(n)得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,包含了窗函数的n个系数。
4.海明窗:
利用w=hamming(n)得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,包含了窗函数的n个系数。
它和汉宁窗的主瓣宽度相同,但是它的旁瓣进一步被压低。
5.布拉克曼窗:
利用w=blackman(n)得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,包含了窗函数的n个系数。
它的主瓣宽度是矩形窗主瓣宽度的3倍,为12*pi/N,但是它的最大旁瓣值比主瓣值低57dB。
6.切比雪夫窗:
它是等波纹的,利用函数w=(N,R)方式设计出N阶的切比雪夫2窗函数,函数的主瓣值比旁瓣值高RdB,且旁瓣是等波纹的。
7.巴特里特窗:
利用w=bartlett(n)的形式得到窗函数,其中n为窗函数的长度,而返回值w为一个n阶的向量,包含了窗函数的n个系数。
8.凯塞窗:
利用w=kaiser(n,beta)的形式得到窗函数。
【阅读文献】
【发现问题】
【问题探究】
【仿真程序】
1.矩形窗(RectangleWindow)调用格式:
w=boxcar(n),根据长度n产生一个矩形窗w。
N=30;
K=300;
t=0:
K-1;
w=boxcar(N);
Y=abs(fftshift(fft(w,K)));
plot(t,Y);
2.三角窗(TriangularWindow)调用格式:
w=triang(n),根据长度n产生一个三角窗w。
N=30;
K=300;
t=0:
K-1;
w=triang(N);
Y=abs(fftshift(fft(w,K)));
plot(t,Y);
3.汉宁窗(HanningWindow)调用格式:
w=hanning(n),根据长度n产生一个汉宁窗w。
N=30;
K=300;
t=0:
K-1;
w=hanning(N);
Y=abs(fftshift(fft(w,K)));
plot(t,Y);
4.海明窗(HammingWindow)调用格式:
w=hamming(n),根据长度n产生一个海明窗w。
N=30;
K=300;
t=0:
K-1;
w=hamming(N);
Y=abs(fftshift(fft(w,K)));
plot(t,Y);
5.布拉克曼窗(BlackmanWindow)调用格式:
w=blackman(n),根据长度n产生一个布拉克曼窗w。
N=30;
K=300;
t=0:
K-1;
w=blackman(N);
Y=abs(fftshift(fft(w,K)));
plot(t,Y);
6.恺撒窗(KaiserWindow)调用格式:
w=kaiser(n,beta),根据长度n和影响窗函数旁瓣的β参数产生一个恺撒窗w。
N=30;
K=300;
t=0:
K-1;
w=kaiser(N,50);
Y=abs(fftshift(fft(w,K)));
plot(t,Y);
【研讨题目】
3.利用DFT分析x(t)=e-|t|u(t)的频谱。
【题目分析】首先利用CTFT公式计算其模拟频谱的理论值;然后对其进行等间隔理想采样,得到序列,利用DTFT公式计算采样序列的数字连续频谱理论值,通过与模拟频谱的理论值对比,理解混叠误差形成的原因及减小误差的措施.接下来是对序列进行加窗处理,得到有限长加窗序列,再次利用DTFT公式计算加窗后序列的数字连续频谱,并与加窗前的数字连续频谱进行对比,理解截断误差形成的原因及减小误差的措施;最后是对加窗序列进行DFT运算,得到加窗后序列的DFT值,它是对数字连续频谱进行等间隔采样的采样值.
【仿真结果】
该程序过程清晰、容易理解,程序运行结果如图2所示。
图2(a)是信号的时域波形,图2(b)是对应的幅度谱图。
由于在归一化频率为0.5的地方还有较大幅度,所以对信号进行理想采样后存在较大的混叠失真,表现在图2(d)中归一化频率为-0.5~0.5范围内的波形与图2(b)的波形明显不同。
图2(e)是理想采样序列加矩形窗得到的图形,对应的幅度谱线如图2(f)中实线所示。
该谱线与图2(d)相比有明显的不同(如波动现象),这是加窗截断的结果。
窗口长度越短截断效应会越明显。
对加窗序列进行DFT运算,只要DFT点数大于等于窗口长度,算出的DFT值就是对加窗序列的连续频谱在一个周期内进行的等间隔采样的采样值。
在图2(f)中表现为DFT幅值在加窗序列连续幅谱的谱线上,是连续频谱曲线上的有限个数据点,即所谓栅栏效应。
图2(f)中用虚线画出了连续信号的幅谱理论值,与DFT的幅值对比,存在一定误差,但只要采样频率足够高,时窗长度足够长,FFT点数足够大,得到的DFT值越逼近实际频谱。
【自主学习内容】
【阅读文献】
【发现问题】
【问题探究】
【仿真程序】
clc;closeall;clear;
%CTFT程序,以x(t)=exp(-t)t>=0为例
%利用数值运算计算并绘制连续信号波形
L=4,%定义信号波形显示时间长度
fs=4,T=1/fs;%定义采样频率和采样周期
t_num=linspace(0,L,100);%取若干时点,点数决定作图精度
xt_num=exp(-1*t_num);%计算信号在各时点的数值
subplot(3,2,1);plot(t_num,xt_num),%绘信号波形
xlabel('时间(秒)'),ylabel('x(t)'),%加标签
grid,title('(a)信号时域波形'),%加网格和标题
%利用符号运算和数值运算计算连续信号幅度谱的理论值
symstW%定义时间和角频率符号对象
xt=exp(-1*t)*heaviside(t),%连续信号解析式
XW=fourier(xt,t,W),%用完整调用格式计算其傅氏变换
%在0两边取若干归一化频点,点数决定作图精度
w1=[linspace(-0.5,0,50),linspace(0,1.5,150)];
XW_num=subs(XW,W,w1*2*pi*fs);%利用置换函数求频谱数值解
mag1=abs(XW_num);%计算各频点频谱的幅值
subplot(3,2,2);plot(w1,mag1),%绘制归一化频率幅值谱线
xlabel('频率(*2*pi*fs)rad/s'),ylabel('幅度'),%加标签
grid,title('(b)连续信号幅频理论值'),%加网格和标题
%DTFT程序,以x(n)=exp(-nT)n>=0为例
%利用数值运算计算并绘制离散信号图形
N=L*fs+1;n_num=0:
N-1;%生成信号波形采样点
xn_num=exp(-1*n_num*T);%计算信号理想采样后的序列值
subplot(3,2,3);stem(n_num,xn_num,'b.'),%绘序列图形
xlabel('n'),ylabel('x(n)'),%加标签
grid,title('(c)理想采样图形'),%加网格和标题
%利用符号运算和数值运算计算离散信号幅度谱的理论值
symsnzw%定义符号对象
xn=exp(-n*T),%定义离散信号
Xz=ztrans(xn,n,z),%用完整调用格式计算其Z变换
%利用复合函数计算序列傅里叶变换的解析解
Z=exp(j*w);Hejw=compose(Xz,Z,z,w);
Hejw_num=subs(Hejw,w,w1*2*pi);%求频谱数值解
mag2=abs(Hejw_num);%计算各频点频谱的幅度
subplot(3,2,4);plot(w1,mag2*T),%绘制频谱幅度曲线
xlabel('频率(*2*pi)rad'),ylabel('幅度'),%加标签
grid,title('(d)离散信号幅频理论值'),%加网格和标题
%序列加窗图示及频谱幅值绘制
%利用数值运算计算并绘制加窗后序列xw(n)的图形
M=8;win=(window(@rectwin,M))';%定义窗点和窗型
xwn=xn_num.*[win,zeros(1,N-M)];%给离散信号加窗
subplot(3,2,5);stem(n_num,xwn,'b.'),%加窗序列图示
xlabel('n'),ylabel('xw(n)'),%加标签
grid,title('(e)加窗序列图形'),%加网格和标题
%利用符号运算和数值运算计算加窗序列的频谱幅值
%先求加窗序列的Z变换,注意表达式长度限制问题
Xwz=0;forn=0:
(M-1);Xwz=Xwz+xwn(n+1)*z^(-n);end
%利用复合函数计算加窗序列傅里叶变换的解析解
Zw=exp(j*w);HejwM=compose(Xwz,Zw,z,w);
HejwM_num=subs(HejwM,w,w1*2*pi);%求频谱数值解
mag3=abs(HejwM_num);%计算各频点频谱的幅度
subplot(3,2,6);plot(w1,mag3*T),%绘频谱幅度曲线
%利用DFT计算加窗序列xw(n)的离散谱幅值
Ndft=16,Xk=fft(xwn,Ndft);%定义DFT点数和DFT运算
Xk0=fftshift(Xk)*T;%将DFT值0对称和幅值加权处理
ifmod(Ndft,2)==0;N1=Ndft;elseN1=Ndft-1;end;
k=[0:
(Ndft-1)]-N1/2;wk=k/Ndft;%0对称取值并归一化
holdon;stem(wk,abs(Xk0),'r.'),%绘制DFT图形
legend('幅谱','DFT',0),%加响应图例,位置自动最佳
xlabel('归一化频率'),ylabel('幅度'),%加标签
grid,title('(f)加窗序列幅谱及其DFT幅值'),
plot(w1,mag1,'k:
'),%与连续信号幅谱的理论值比较
【研讨题目】
4.利用DFT分析计算x(t)=cos2pit的频谱,设w0=pi/2
【题目分析】因为f=w/2pi=1Hz,所以抽样频率需要大于2Hz.
【序列频谱计算的基本方法】
【仿真结果】
【结果分析】
【自主学习内容】已知一连续信号为
其中
。
试利用DFT分析其频谱。
解:
信号
的最高频率
,抽样频率
,取抽样频率
;最低的频率分辨率为
,最少的信号样点数为
。
的MATLAB程序如下:
%programexa_1_1.m,利用矩形窗计算有限长余弦信号频谱
N=30;%数据的长度
L=512;%DFT的点数
f1=100;
f2=120;
fs=600;%抽样频率
T=1/fs;%抽样间隔
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
y=fft(x,L);
mag=abs(y);
f=(0:
length(y)-1)'*fs/length(y);
plot(f(1:
L/2),mag(1:
L/2));
xlabel('频率(Hz)')
ylabel('幅度谱')
程序运行结果如下图所示。
由图可见,频谱图显示出两个较为明显的峰值(对应
)。
结论:
当截取信号
样点时,频率分辨率
,刚好能够分辨出
和
两个频谱分量,但频谱泄漏较严重。
若取
(达不到最低的频率分辨率)。
的MATLAB程序如下:
%programexa_1_2.m,利用矩形窗计算有限长余弦信号频谱
N=20;%数据的长度
L=512;%DFT的点数
f1=100;
f2=120;
fs=600;%抽样频率
T=1/fs;%抽样间隔
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
y=fft(x,L);
mag=abs(y);
f=(0:
length(y)-1)'*fs/length(y);
plot(f(1:
L/2),mag(1:
L/2));
xlabel('频率(Hz)')
ylabel('幅度谱')
程序运行结果如下图所示。
由图可见,频谱图显示不出对应
的两个明显峰值。
结论:
当截取信号
样点时,频率分辨率
,达不到最低的频率分辨率
,故分辨不出
和
两个频谱分量,且频谱泄漏更为严重。
若取频率分辨率
,则对应的信号样点数为
。
的MATLAB程序如下:
%programexa_1_3.m,利用矩形窗计算有限长余弦信号频谱
N=600;%数据的长度
f1=100;
f2=120;
fs=600;%抽样频率
T=1/fs;%抽样间隔
t=(0:
N-1)*T;
x=cos(2*pi*f1*t)+cos(2*pi*f2*t);
y=fft(x);
mag=abs(y);
f=(0:
length(y)-1)'*fs/length(y);
plot(f(1:
length(y)/2),mag(1:
length(y)/2));
xlabel('频率(Hz)')
ylabel('幅度谱')
程序运行结果如下图所示。
由图可见,频谱图显示出两个特别明显的峰值(对应
)。
结论:
当截取信号
样点时,频率分辨率
,高分辨率的频谱图具有较高的质量,频谱分析时必须保证获取足够的信号数据长度。
【阅读文献】
【发现问题】
【问题探究】
【仿真程序】
clc;closeall;clear;
%CTFT程序,以x(t)=exp(-t)t>=0为例
%利用数值运算计算并绘制连续信号波形
L=8,%定义信号波形显示时间长度
fs=4,T=1/fs;%定义采样频率和采样周期
t_num=linspace(0,L,100);%取若干时点,点数决定作图精度
xt_num=cos(2*pi*t_num);%计算信号在各时点的数值
subplot(3,2,1);plot(t_num,xt_num),%绘信号波形
xlabel('时间(秒)'),ylabel('x(t)'),%加标签
grid,title('(a)信号时域波形'),%加网格和标题
%利用符号运算和数值运算计算连续信号幅度谱的理论值
symstW%定义时间和角频率符号对象
xt=cos(2*pi*t)*heaviside(t),%连续信号解析式
XW=fourier(xt,t,W),%用完整调用格式计算其傅氏变换
%在0两边取若干归一化频点,点数决定作图精度
w1=[linspace(-0.5,0,50),linspace(0,1.5,150)];
XW_num=subs(XW,W,w1*2*pi*fs);%利用置换函数求频谱数值解
mag1=abs(XW_num);%计算各频点频谱的幅值
subplot(3,2,2);plot(w1,mag1),%绘制归一化频率幅值谱线
xlabel('频率(*2*pi*fs)rad/s'),ylabel('幅度'),%加标签
grid,title('(b)连续信号幅频理论值'),%加网格和标题
%利用数值运算计算并绘制离散信号图形
N=L*fs+1;n_num=0:
N-1;%生成信号波形采样点
xn_num=cos(2*pi*n_num*T);%
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- DSP 研究性学习 报告