1、实验名称：频谱峰值搜索实验目的：1、熟练掌握离散信号的 DFT 实现方法；2、熟练掌握 Matlab 实现 DFT 的方法，提高编程实践能力；3、增强自我学习能力，查阅文献搜索能力；4、掌握离散信号的时域与频域的对应关系。实验原理：1. 离散复正弦信号的 DFT(1)210()()NjnknXkfxe2、Matlab 主要函数fft(signal,N); signal:输入信号， N：fft 的点数 函数的作用是对输入信号做 N 点的 DFTfftshift(fft(signal,N); 将零频点移到频谱的中间3 谱峰搜索算法采用一维黄金分割精搜算法5取代分级搜索过程中的递归精搜。对于一维黄金
2、分割精搜算法, 若函数 有且仅有一个极大值位于区间 上, 则有当()fx()ac时, 且 。此时若在区间 上选取一点 , 当abc()fabfc,bx时 , 则 且 , 即极大值点将位于三元点组()fx()f()x对应的区间 上； 否则, 当 时, 则 且,x()f()fb, 即极大值点将位于三元点组 对应的区间 上。()()fbfcbc,c在这些三元点组中, 其中间点对应的函数值都是每一轮求解过程中的最大值.这一过程下图所示, 继续对三元点组所对应的区间进行划分 , 直到区间足够小, 小到以至于该区间上任何一点都可以表示函数的极大值点。下面是一维黄金分割搜索算法的基本原理。给定三元点组 ,
3、假设 是 , 之间的一个分割 。即abcacw(2)1bw再假设一个试探点 位于区间 之间, 且有x()bc（3） xzca通过选取试探点 后, 可以将极大值点压缩到相对长度为 的区间 。或x wz(,)ax者是相对长度为 的区间 上。考虑到搜索极大值的最坏情况,应该使得1w()b下式成立:（4）1 z=-2zw可见试探点 应该选为点 关于区间 的对称点。同样应该保证 是xb(,)acx之间的一个分割 。(,)bcw（5）1xzcbw把式( 4) 代入到式 (5) 中得到如下的二次方程（6）230解得其根为( 考虑到 , 舍弃另一根)01w（7）5.81962所以对于给定的三元点组所对应的区间
4、上, 每次选取试探点都是位于较大的一段子区间上, 并且距离原来中间点0.3189660 的位置上。也就是说通过一次这样的试探点的选取可以使得极大值点将位于原区间0.6180340 的更小的区间上。继续这样的计算, 那么区间将变得越来越小, 小到以至于该区间上任何一点都可以表示函数的极大值点。实验步骤：1、设置输入信号的参数以及 DFT 变换的点数根据要求，输入信号的模拟频率为 ，那么采样频率满足10.1f即可，为方便观察频率最大值位置，取 =2Hz。给定 DFT 点数为 102412sf s点，而为了使的被观察的频谱峰值在频谱图的中央，将抽样时间取在 1,ssNtff的区间，采样间隔为 。其中
5、 N=512，满足采样点数为 点。1/ssTf 51204这样得到输入信号的表达式为（8）12sjfnTigale2、对信号进行 DFT，并画出频谱图。（1）在 MATLAB 中应用 fft（signal，N）对信号 signal 做 N 点的 FFT；（2）分别应用函数 fftshift、abs 对 DFT 结果调整和取绝对值；（3）设置横坐标。根据 和 的对应关系，得到很坐标的取值范围是-sf1,1。在 MATLAB中设置为f=(1:2*N)-N)*(fs/(2*N);3、运用一维黄金分割方法找出频谱峰值。根据频谱的峰值范围，以及分割法的原理，设置个参数为： ，0.5a， 。创建计算相应频
6、率点的幅值计算函数。根据 DFT 的计0.5c.31896w算定义， 处的频谱值为xf（9）221100()()()xsfNNjnNjnkxssnfTexTe根据计算精度，将MATLAB计算精度设置为format long 。并设计计算迭代次数的变量iterations。通过判别 b点和x点的幅值大小来更新参数，参数更新如下（10）() ,abxcf则 ：则 ：实验结果：1、运行程序（程序见附录） ，得到频谱图如图 1 所示图 1 复正弦信号的频谱由于图上显示精度的原因，直接找到的最大值不是我们所需要的最大值，通过峰值搜索函数得到最大值。2、得到搜索结果为：iterations = 53。ma
7、x0.147039f实验结果分析：1、由于 DFT 的点数 1024 比较多，而频谱范围较小，所以离散的频谱在图上显得像连续谱一样。2、因为输入信号的模拟频率的值为 =0.111111111，所以图上离散的点上没有1f显示最大值点。3、从搜索的结果看出，最大值在满足精度要求的情况下是正确的，说明一维黄金搜索方法在本实验中是可行的。4、在 取值距离最大值较远的情况下经过 53 次循环迭代可以得到最大值。,ac说明该算法收敛比较快。附录：clear all;clc;format long N=512;设置采样点数为2N=1024fs=2;%设置采样频率为2Hzt=-N/fs:1/fs:N/fs-1
8、/fs;%采样时间序列f1=0.111111111;输入信号的频率signal=exp(i*2*pi*f1*t);输入信号的采样序列signalDtf=abs(fftshift(fft(signal); %对信号进行DFT%画出波形f=(1:2*N)-N)*(fs/(2*N);设置横坐标plot(f,20*log10(signalDtf);hold on;gridxlabel(f Hz);ylabel(20log10(幅度);title(输入信号的DFT);%采用一维黄金分割精度算法%结合算法的特点选择a=-0.2,c=-0.2,w=0.3189660;%初始化幅值ampX=2;ampB=0;
9、a=-0.5;c=0.5;w=0.3189660;z=1-2*w;iterations=0;%初始化迭代次数while abs(ampX-ampB)1.0000e-0011;%设置收敛目标%定义算法参数 b=w*(c-a)+a;x=z*(c-a)+b;iterations=iterations+1;迭代次数更新%通过DFT 定义计算a，b，c ，x 点的幅值ampB=Amplitude(b);ampA=Amplitude(a);ampC=Amplitude(c);ampX=Amplitude(x);compairBX=ampBampX;%比较b，x 点的幅值大小if compairBX=1;a=a;b=b;c=x;else compairBX=0;a=b;b=x;c=c;endendamplitudeMax=xiterationsfunction ampx=Amplitude(x)%频率幅值对应幅值计算函数%计算任意频率对应点的幅度值%变量x为模拟平频率%以fs=1 的采样频率对其进行采样% x 频率值format long N=1024;fs=1;n=0:N-1;f1=0.111111111;signal=exp(j*2*pi*f1*n);wight=exp(j*2*pi/N)*n*x*N/fs);ampx=20*log10(abs(signal*wight);end
