频谱峰值搜索

实验名称： 频谱峰值搜索
实验目的：
1、熟练掌握离散信号的DFT 实现方法；
2、熟练掌握Matlab 实现DFT 的方法，提高编程实践能力；
3、增强自我学习能力，查阅文献搜索能力；
4、掌握离散信号的时域与频域的对应关系。

实验原理：
1. 离散复正弦信号的DFT
21
10()()N j nk N
n X kf x n e
π-==∑ (1)
2、Matlab 主要函数
fft(signal,N); signal:输入信号，N ：fft 的点数 函数的作用是对输入信号做N 点的DFT fftshift(fft(signal,N)); 将零频点移到频谱的中间 3谱峰搜索算法
采用一维黄金分割精搜算法[5]取代分级搜索过程中的递归精搜。

对于一维黄金分割精搜算法, 若函数()f x 有且仅有一个极大值位于区间(,)a c 上, 则有当
a b c <<时,()()f b f a > 且()()f b f c >。

此时若在区间(,)b c 上选取一点x , 当
()()f b f x >时, 则()()f b f a >且()()f b f x > , 即极大值点将位于三元点组
a b x <<对应的区间(,)a x 上； 否则, 当()()f b f x < 时, 则()()f x f b >且
()()()f x f b f c >>, 即极大值点将位于三元点组b x c << 对应的区间(,)b c 上。

在这些三元点组中, 其中间点对应的函数值都是每一轮求解过程中的最大值.这一过程下图所示, 继续对三元点组所对应的区间进行划分, 直到区间足够小, 小到以至于该区间上任何一点都可以表示函数的极大值点。

下面是一维黄金分割搜索算法的基本原理。

给定三元点组a b c <<, 假设b 是a ,c 之间的一个分割w 。

即
, 1b a c b
w w c a c a --==--- (2)
再假设一个试探点x 位于区间(,)b c 之间, 且有
x b
z c a -=- （3）
通过选取试探点x 后, 可以将极大值点压缩到相对长度为w z +的区间(,)a x 。

或者是相对长度为1w -的区间(,)b c 上。

考虑到搜索极大值的最坏情况,应该使得下式成立:
1 z=1-2w w z w +=- （4）
可见试探点x 应该选为点b 关于区间(,)a c 的对称点。

同样应该保证x 是
(,)b c 之间的一个分割w 。

1x b z
w w c b w -==-- （5）
把式( 4) 代入到式(5) 中得到如下的二次方程
2310w w -+= （6）
解得其根为( 考虑到01w <<, 舍弃另一根)
30.38196602
w =
= （7） 所以对于给定的三元点组所对应的区间上, 每次选取试探点都是位于较大的一段子区间上, 并且距离原来中间点0.3189660 的位置上。

也就是说通过一次这样的试探点的选取可以使得极大值点将位于原区间0.6180340 的更小的区间上。

继续这样的计算, 那么区间将变得越来越小, 小到以至于该区间上任何一点都可以表示函数的极大值点。

实验步骤：
1、设置输入信号的参数以及DFT 变换的点数
根据要求，输入信号的模拟频率为10.111111111f =，那么采样频率满足
12s f f >即可，为方便观察频率最大值位置，取s f =2Hz 。

给定DFT 点数为1024
点，而为了使的被观察的频谱峰值在频谱图的中央，将抽样时间取在
1[
,]s s s
N N t f f f =- 的区间，采样间隔为1/s s T f =。

其中N=512，满足采样点数为51221024⨯=点。

这样得到输入信号的表达式为
12s j f nT signal e π= （8）
2、对信号进行DFT ，并画出频谱图。

（1）在MATLAB 中应用fft （signal ，N ）对信号signal 做N 点的FFT ； （2）分别应用函数fftshift 、abs 对DFT 结果调整和取绝对值；
（3）设置横坐标。

根据s f 和N 的对应关系，得到很坐标的取值范围是[-1,1]。

在MATLAB 中设置为f=((1:2*N)-N)*(fs/(2*N)); 3、运用一维黄金分割方法找出频谱峰值。

根据频谱的峰值范围，以及分割法的原理，设置个参数为：0.5a =-，0.5c =，0.318966w =。

创建计算相应频率点的幅值计算函数。

根据DFT 的计算定义，x f 处的频谱值为
221
1
()()()x s
f N N j
nN j
nk N f N
x s s n n x f x nT e
x nT e
π
π--====∑∑ （9）
根据计算精度，将MATLAB 计算精度设置为format long 。

并设计计算迭代次数的变量iterations 。

通过判别b 点和x 点的幅值大小来更新参数，参数更新如下
()() ,,()() ,,x b x x a a b b x c
x b x f b a x b c c
>→→→⎧⎨
<→→→⎩则：则： （10）
实验结果：
1、运行程序（程序见附录），得到频谱图如图1所示
图 1 复正弦信号的频谱
由于图上显示精度的原因，直接找到的最大值不是我们所需要的最大值，通过峰值搜索函数得到最大值。

2、得到搜索结果为：
max 0.111111111474039f = iterations = 53。

实验结果分析：
1、由于DFT的点数1024比较多，而频谱范围较小，所以离散的频谱在图上显得
像连续谱一样。

2、因为输入信号的模拟频率的值为
f=0.111111111，所以图上离散的点上没有
1
显示最大值点。

3、从搜索的结果看出，最大值在满足精度要求的情况下是正确的，说明一维黄
金搜索方法在本实验中是可行的。

4、在,a c取值距离最大值较远的情况下经过53次循环迭代可以得到最大值。

说
明该算法收敛比较快。

附录：
clear all;clc;format long
N=512;设置采样点数为2N=1024
fs=2;%设置采样频率为2Hz
t=-N/fs:1/fs:N/fs-1/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;grid
xlabel('f Hz');ylabel('20log10(幅度)');title('输入信号的DFT');
%采用一维黄金分割精度算法
%结合算法的特点选择a=-0.2,c=-0.2,w=0.3189660;
%初始化幅值
ampX=2;ampB=0;
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=ampB>ampX;%比较b，x点的幅值大小
if compairBX==1;
a=a;b=b;c=x;
else compairBX==0;
a=b;b=x;c=c;
end
end
amplitudeMax=x
iterations
function 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。