【阅读文献】
[1]陈后金,薛健,胡健.数字信号处理[M].北京:
高等教育出版社,2006.
【发现问题】(专题研讨或相关知识点学习中发现的问题):
fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。
因为直接用fft得出的数据与频率不是对应的,fftshift可以纠正过来
【问题探究】
在离散序列频谱计算中为何要用窗函数?
用不同的窗函数对计算结果有何影响?
与矩形窗相比哈明窗有何特点?
如何选择窗函数?
答:
有一些离散序列长度无限长,计算机无法处理,所以要利用窗函数进行截短。
用不同的窗函数得到的计算结果不同,与矩形窗相比哈明窗减小了旁瓣,却加宽了主瓣宽度。
【仿真程序】
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);
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');
3已知信号
,
(1)以
为间隔对
抽样,抽得128个样本
,画出其频谱密度函数的草图
,对其做128点FFT,得到
,问m为哪些值时,
具有较大值?
(较大值即主瓣之中的非零值)
(2)若给
补上896个零值,使之成为1024点序列,
对
做1024点FFT得到
,问m为哪些值时,
具有较大值?
(3)仍以
对
抽样,抽得1024个样本
,画出其频谱密度函数的草图
,对
做1024点FFT,得到
问m为哪些值时,具
有较大值?
(4)设
为1024点的hamming窗,
写出
的数学表达式,令
,画出
的频谱密度函数的草图
,对
做1024点FFT,得到
,问m为哪些值时,
具有较大的值?
【题目分析】
分析影响谱峰分辨率的主要因数,进一步认识补零在在频谱计算中的作用。
【仿真结果】
(1)
从上图看出,m=51、52、53及75、76、77时,
具有较大值
(2)
从上图看出,m=403—415及609—621时,
具有较大值
(3)
从上图看出,m=403—415及609—621时,
具有较大值
(4)
从上图看出,m=403—415及609—621时,
具有较大值
【结果分析】
DFT点数越多,分辨率越高,补零只能使序列的频谱变得细致,但不能提高频率的分辨率。
只有采集更多的有效数据,才能得到序列的高分辨频谱。
【自主学习内容】
【阅读文献】
[1]陈后金,薛健,胡健.数字信号处理[M].北京:
高等教育出版社,2006.
【发现问题】(专题研讨或相关知识点学习中发现的问题):
【问题探究】
1、2、3题讨论的是离散信号频谱的计算问题。
与连续信号频谱计算问题相比较,其计算误差有何不同?
【仿真程序】
(1)
f1=10*10^3;
f2=10.3*10^3;
T=0.01*10^-3;
L=128;
k=0:
L-1;
t=T*k;
x1=cos(2*pi*f1*t)+cos(2*pi*f2*t);
X1=fftshift(fft(x1));
plot((-L/2:
L/2-1)*2*pi/L,abs(X1))
xlabel('\Omega')
ylabel('X1(exp(j\Omega))')
title([num2str(L)'点FFT'])
[ab]=max(X1)
(2)
f1=10*10^3;
f2=10.3*10^3;
T=0.01*10^-3;
N=128;
k=0:
N-1;
L=1024;
t=T*k;
x1=cos(2*pi*f1*t)+cos(2*pi*f2*t);
G1=fftshift(fft(x1,L));
stem(-L/2:
L/2-1,abs(G1))
xlabel('k')
ylabel('|G1[m]|')
title([num2str(L)'点FFT'])
[ab]=max(G1)
(3)
f1=10*10^3;
f2=10.3*10^3;
T=0.01*10^-3;
N=1024;
k=0:
N-1;
L=1024;
t=T*k;
x2=cos(2*pi*f1*t)+cos(2*pi*f2*t);
X2=fftshift(fft(x2,L));
plot((-L/2:
L/2-1)*2*pi/L,abs(X2))
xlabel('\Omega')
ylabel('X2(exp(j\Omega))')
title([num2str(L)'点FFT'])
(4)
f1=10*10^3;
f2=10.3*10^3;
T=0.01*10^-3;
N=1024;
k=0:
N-1;
L=1024;
t=T*k;
x2=cos(2*pi*f1*t)+cos(2*pi*f2*t);
wh=(hamming(N))';
x2=x2.*wh;
G2=fftshift(fft(x2,L));
plot((-L/2:
L/2-1)*2*pi/L,abs(G2))
xlabel('\Omega')
ylabel('G2(exp(j\Omega))')
title([num2str(L)'点FFT'])
4试用DFT近似计算高斯信号
的频谱抽样值。
高斯信号频谱的理论值为
通过与理论值比较,讨论信号的时域截取长度和抽样频率对计算误差的影响。
【题目分析】
连续非周期信号频谱计算的基本方法。
计算中出现误差的主要原因及减小误差的方法。
对于连续非周期信号,要对其进行频谱计算,需要先经过抽样,变为离散信号,再对离散信号的频谱在一个周期内抽样,即做DFT。
由于题目所给信号时域无线长,所以为了在抽样后频谱不混叠,对信号进行截短即加抗混叠滤波器后在进行频谱分析,因为加滤波器后照成的误差远小于频谱混叠照成的误差。
在对信号抽样时还需要考虑抽样频率,以使抽样后所得频谱不会发生混叠。
【仿真结果】
【结果分析】
由以上三幅图可以看出:
当时域截取长度相同时,抽样间隔越小时误差越小,当抽样间隔一定时,时域截取长度越长,误差越小。
当取抽样间隔为1S,时域截取长度为2S时,误差较大,绝对误差在0.5左右;当抽样间隔为0,5S,时域截取长度为2S时,误差比间隔为1S时小,绝对误差不大于0.2;当抽样间隔为0.5S时域截取长度为4S时,误差更小,绝对误差不大于0.04。
因为时域截取长度越长,保留下来的原信号中的信息越多,抽样间隔越小,频谱越不容易发生混叠,所以所得频谱与理论值相比,误差更小。
【自主学习内容】
【阅读文献】
[1]陈后金,薛健,胡健.数字信号处理[M].北京:
高等教育出版社,2006.
【发现问题】(专题研讨或相关知识点学习中发现的问题):
【问题探究】
经过对比,当我们取值Ts=0.05,N=40时,时域截取长度为2s,其产生的实验误差绝对值小于0.01;当我们Ts=0.05,N=64时,时域截取长度为3.2s,实验误差绝对值被控制在0.0005以内。
但由上图可以看到,Ts=0.05,N=100时的实验误差绝对值甚至大于0.2.
【仿真程序】
Ts=0.5;
N=4;
N0=64;
k=(-N/2:
(N/2))*Ts;
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(j\omega)');
legend('试验值','理论值');
title(['Ts=',num2str(Ts)'''N=',num2str(N)]);
subplot(2,1,2)
plot(w/pi,abs(X)-XT)
ylabel('实验误差')
xlabel('\omega/\pi');
扩展题
5本题研究连续周期信号频谱的近似计算问题。
周期为T0的连续时间周期信号x(t)可用Fourier级数表示为
其中
X(nω0)称为连续时间周期信号x(t)的频谱函数。
称为信号的基频(基波),
称为信号的谐波。
如果信号x(t)函数表达式已知,则可由积分得出信号的频谱。
如果信号x(t)函数表达式未知,或者x(t)函数表达式非常复杂,则很难由积分得信号的频谱。
本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。
(1)若在信号x(t)的一个周期T0内抽样N个点,即
,T为抽样周期(间隔),可获得序列x[k]
试分析序列x[k]的DFT与连续时间周期信号x(t)的频谱X(nω0)的关系;
(2)由
(1)的结论,给出由DFT近似计算周期信号频谱X(nω0)的方案;
(3)周期信号x(t)的周期T0=1,x(t)在区间[0,1]的表达式为
x(t)=20t2(1-t)4cos(12πt)
(a)试画出信号x(t)在区间[0,1]的波形;
(b)若要用6次以内的谐波近似表示x(t),试给出计算方案,并计算出近似表示的误差。
讨论出现误差的原因及减小误差的方法。
【题目分析】
【理论推导】
DFT计算所得结果X[m]与连续周期信号频谱X(nω0)的关系。
DFT计算所得结果X[m]与连续周期信号频谱X(nω0)的关系。
由于
于是得
令n=m+rN;m=0,1,2,…N-1,r为整数,上式可化为
化简可得
对比IDFT
可得
【计算方案】
根据理论推导结果设计近似计算方案。
分析产生误差的主要原因。
从上面的关系式可以看出X[m]是连续轴系信号频谱X(now)发生周期化,再乘上比例因子N所得,由于对于一般的连续周期信号的频谱,当n→±∞时,┃X(nwo0┃→0,所以用X[m]来近似为较低的谐波分量,当然这一过程中的误差很显然主要在于低次谐波分量与高次谐波分量发生的混叠。
如果在一个周期内取较多的抽样点N,可以提高与低频分量发生混叠的高频分量的频率值,高频分量如果足够高,也就是幅度值可以趋近于零,就可以进行有效的近似,其混叠误差也会有效地减少。
同时,在进行信号的重建时,谐波分量选取的次数的多少对误差也有影响。
对于一般的周期信号,通常时域是无限的,而我们只能选择次数较低的谐波分量(占有原信号的主要能量)进行信号重建,这是不可避免的会发生泄漏,产生误差。
【扩展分析】
如果周期信号x(t)是带限信号,即信号的最高频率分量为Mω0(是正整数),试确定在一个周期内的最少抽样点N,使得在频谱的计算过程当中不存在混叠误差。
与抽样定理给出的结论比较,发表你的看法。
由推出的理论关系式,频谱发生的是Now的周期化,原信号x(t)带限,最高频率分量为Mwo,考虑到想x(t)是实信号,一定也会有-Mwo分量,当N=2M时,频谱会发生混叠,所以选取N>2M,这是再将频谱周期化将不会发生混叠,这与抽样定理相似。
这个结果是抽样定理的结果,只是周期信号是一种特殊信号,当然也符合抽样定理。
只是不能取临界条件。
【仿真结果】
一个周期内抽样点个数不同,结果比较:
N=32时最大误差在
左右,平均误差不到
,而N=16时,最大误差接近0.15,平均误差接近0.07,可以看出N=32时的误差比N=16时的误差小得多,说明N=16时发生的混叠严重。
选取谐波分量的次数对实验近似的影响:
选取8次谐波分量,误差在0.03左右,平均误差在0.015左右,
可以看出当只取8次谐波分量重建新号,已经有了相当的误差,只取六次谐波分量重建新号,误差已经很大了,所以选取合适次数的谐波进行信号重建也对信号是否能合理恢复有着很大的影响。
【结果分析】
讨论DFT点数对近似计算的影响,讨论所取谐波项的多少对近似计算的影响。
误差分析要给出定量的结果,如平均误差,最大误差等。
与连续非周期信号频谱计算过程中存在的误差相比较,连续周期信号频谱的计算计算误差有何异同?
【自主学习内容】
【阅读文献】
[1]陈后金,薛健,胡健.数字信号处理[M].北京:
高等教育出版社,2006.
【发现问题】(专题研讨或相关知识点学习中发现的问题):
d值的选取问题。
由理论计算公式可知,d值选为pi更为方便。
同时,公式可知,进行DFT变换之前函数还要乘T,这样才能与G(jw)对应。
【问题探究】
【仿真程序】
T0=1;
N=32;
T=T0/N;
n=11;
ts=(0:
N-1)*T;
x=20*ts.^2.*((1-ts).^4).*cos(12*pi.*ts);
X=fft(x);
c=[X(N-n+2:
N)X(1:
n)]/N;
t=linspace(0,1,1000);
xt=20*t.^2.*((1-t).^4).*cos(12*pi.*t);
xe=zeros(1,length(t));
fork=-(n-1):
(n-1),
xe=xe+c(k+n)*exp(j*k*2*pi*t/T0);
end
figure
(1)
plot(t,xt,t,xe);
gridon;
title(['抽取点数N=',num2str(N)'时的近似情况']);
xlabel('t(s)');
ylabel('x(t)');
legend('理论值','近似值');
figure
(2)
plot(t,xt-xe);
gridon;
title(['抽取点数N=',num2str(N)'时的误差情况']);
xlabel('t(s)');
ylabel('x(t)');
legend('误差');