哪位大神帮小弟写一段MATLAB提取语音信号共振峰程序

function [fmt] = seekfmts1(sig,Nt,fs,Nlpc)

if nargin<4, Nlpc = round(fs/1000)+2end

ls=length(sig) % 数据长度

Nwin = floor(ls/Nt) % 帧长

for m=1:Nt,

lpcsig = sig((Nwin*(m-1)+1):min([(Nwin*m) ls]))% 取来一帧信号

if ~isempty(lpcsig),

a = lpc(lpcsig,Nlpc)% 计算LPC系数

const=fs/(2*pi) % 常数

rts=roots(a)% 求根

k=1 % 初始化

yf = []

bandw=[]

for i=1:length(a)-1

re=real(rts(i)) % 取根之实部

im=imag(rts(i)) % 取根之虚部

formn=const*atan2(im,re)% 计算共振峰频率

bw=-2*const*log(abs(rts(i)))% 计算带宽

if formn>150 &bw <700 &formn<fs/圆冲2 % 满足条件方能成共振峰和带宽

yf(k)=formn

bandw(k)=bw

k=k+1

end

end

[y, ind]=sort(yf) % 排序

bw=bandw(ind)

F = [NaN NaN NaN] % 初始化

F(1:min(3,length(y))) = y(1:min(3,length(y)))% 输出最多三个

F = F(:)% 按列输出

fmt(:,m)=F/(fs/2) % 归一化频率

end

end

from 《matlab在语音信号猜腔慎分析与穗敬合成中的应用》

有很多中文字符，修改下:

%%%%%%%%%%%%%%%%%%%%%%%%

clear %消除内存中所有变量和函数

clc %消除工作窗口所显示的内容

close all hidden%关闭所有隐藏的窗口

%%%%%%%%%%%%%%%%%%%%%%%%%

%提示用键盘输入输入数据文件名

fni=input('消除多项式趋势项-输入数据文件名','s')

%以只读方式打开数据文件

fid=fopen(fni,'r')

sf=fscanf(fid,'%f',1) %读入采样频率值

m=fscanf(fid,'%d',1)%读入拟合多项式阶数

fno=fscanf(fid,'态返%s',1) %读入输出数据文件名

x=fscanf(fid,'%f',inf) %读入时程数据存成列向量

%关闭数据文件

status=fclose(fid)

%取信号数据长度

n=length(x)

%建立离散时间列向量

t=(0:1/sf:(n-1)/sf)

%计算趋势项的多项式系数待定系数向量a

a=polyfit(t,x,m)

%用x减去多项式系数a生成的趋势项

y=x-polyval(a,t)

%将分成2行1列的图形窗口的第1列设为当前绘图区域

subplot(2,1,1)

%绘制x对于t的时程曲线图形

plot(t,x)

%在图幅上添加坐标网络

grid on

%将分成2行1列的图形窗口的第2列设为当前绘图区域

subplot(2,1,2)

%绘制y对于t的时程曲线

plot(t,y)

%在图幅上添加坐标网络

grid on

%以写的方式打开文件或建立一个新文件

fid=fopen(fno,'w')

%进行n次循环将计算结果烂仔写到输出数据文件中

for k=1:n

%每行输出两个实型数据帆历饥，t为时间,y为消除趋势项后的结果

fprintf(fid,'%f%f\n',t(k),y(k))

%循环体结束语句

end

%关闭数据文件

status=fclose(fid)

可能还是没修改完全 你再查查看

好眼熟啊。。。通院备孙08级的随机信号实验题目？

这份是我在网备滚卖络上找的：

（1）.周期图法： 思想：周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度 N。由于序列 x(n)的离散傅里叶变换 X(k) 具有周期性，因而这种功率谱也具有周期性，常称为周期图。周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。

程序：

%首先，生成输入信号的程序为：

clear fs=20000 n=0:1/ fs: 0. 1N=lengt h( n)W=2000*pi%因方波频率 F=1000HZ 所以角频率 W=2000pi X1n=square( W*n) %方波信号 X2n=randn( 1, N)%白噪声信号 xn=X1n+0. 2*X2nsubplot (3, 1, 1) plot (n, xn) xlabel( ' n' ) ylabel('输入信仿逗号') %其次，开始用周期图法进行估计； clear all fs=20000n=0:1/fs:0.1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nsubplot(2,1,1) plot(n,xn)Nfft=256N=256%傅里叶变换的采样点数256 Pxx=abs(fft(xn,Nfft).^2)/Nf=(0:length(Pxx)-1)*fs/length(Pxx)subplot(2,1,2), plot(f,10*log10(Pxx)),%转成DB 单位

（2）.自相关函数法： 思想： 随机信号 x(n)的相关函数是在时间域内描述随机过程的重要特征。自相关函数是随机信号在不同时刻的值之间的依赖性的量度，是一个很有用的统计平均量，其定义为自相关函数 (1) 式中E【·】表示数学期望，*表示共轭值,m 为时间滞后数。在随机信号处理中，自相关函数可以用来检测淹没在随机噪声干扰中的信号，随机信号的自功率谱等于它的自相关函数的傅里叶变换。因此，通过自相关估计可求得信号的功率谱。利用计算机计算自相关估值有两种方法。一种是直接方法，先计算出随机信号和它的滞后序列的乘积，再取其平均值即得相关函数的估计值。另一种是间接方法，先用快速变换算法计算随机序列的功率谱密度，再作反变换计算出相关函数。

程序： n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(2) subplot(3,1,1) plot(n,xn)%输入信号 axis([0 0.01 -5 5]) m =-100:100 [r,lag]=xcorr(xn,100,'biased')%求XN 的自相关函数R，biased 为有偏估计,lag 为R 的序列号 subplot(3,1,2) hndl=stem(m,r)%绘制离散图，分布点从-100—+100 set(hndl,'Marker','.') set(hndl,'MarkerSize',2)ylabel('自相关函数R(m)') %利用间接法计算功率谱 k=0:1000%取1000 个点 w=(pi/500)*kM=k/500X=r*(exp(-li*pi/500).^(m'*k))%对R 求傅里叶变换,li 就是j magX=abs(X)subplot(3,1,3) plot(M,10*log10(magX))xlabel('功率谱的改进直接法估计')

（3）.自协方差法： 思想：在实际中,有时用随机信号在两个不同时刻t 1 ,t 2 的取值X(t 1 )和X(t 2 )之间的二阶混合中心矩来描述随机信号 X(t)在任意两个时刻取值起伏变化的相依程度, 即自协方差。自协方差函数与自相关函数描述的特性是一致的，所以其原理与相关函数法近似。

clear allfs=20000n=0:1/fs:0.1P=2000*piy=square(P*n)xn=y+0.2*randn(size(n))%绘制信号波形 figure(3) subplot(2,1,1) plot(n,xn) xlabel('时间（s）') ylabel('幅度') title('y+randn(size(n))') ymax_xn=max(xn)+0.2ymin_xn=min(xn)-0.2axis([0 0.3 ymin_xn ymax_xn]) %使用协方差法估计序列功率谱 p=floor(length(xn)/3)+1nfft=1024[xpsd,f]=pcov(xn,p,nfft,fs,'half')%绘制功率谱估计 pmax=max(xpsd)xpsd=xpsd/pmaxxpsd=10*log10(xpsd+0.000001)subplot(2,1,2) plot(f,xpsd) title('基于协方差的功率谱估计') ylabel('功率谱估计（db）') xlabel('频率') grid onymin=min(xpsd)-2ymax=max(xpsd)+2axis([0 fs/2 ymin ymax])

（4）.最大熵法 思想：上网查的最大熵法原理是取一组时间序列，使其自相关函数与一组已知数据的自相关函数相同，同时使已知自相关函数以外的部分的随机性最强，以所取时间序列的谱作为已知数据的谱估值。它等效于根据使随机过程的熵为最大的原则,利用 N 个已知的自相关函数值来外推其他未知的自相关函数值所得到的功率谱。最大熵法功率谱估值是一种可获得高分辨率的非线性谱估值方法，特别适用于数据长度较短的情况。最大熵法谱估值对未知数据的假定 ，一个平稳的随机序列，可以用周期图法对其功率谱进行估值。这种估值方法隐含着假定未知数据是已知数据的周期性重复。现有的线性谱估计方法是假定未知数据的自相关函数值为零，这种人为假定带来的误差较大。最大熵法是利用已知的自相关函数值来外推未知的自相关函数值，去除了对未知数据的人为假定，从而使谱估计的结果更为合理。熵在信息论中是信息的度量，事件越不确定，其信息量越大，熵也越大。对于上述问题来说，对随机过程的未知的自相关函数值，除了从已知的自相关函数值得到有关它的信息以外，没有其他的先验知识。因而，在外推时，不希望加以其他任何新的限制，亦即使之“最不确定”。换言之，就是使随机过程的熵最大。

程序： fs=20000n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(5) subplot(3,1,1) plot(n,xn)asis([0 0.01 -5 5]) Nfft=256%分段长度256 [Pxx,f]=pmem(xn,14,Nfft,fs)%调用最大熵函数pmem，滤波器阶数14 subplot(2,1,2), plot(f,10*log10(Pxx)), title(' 最大熵法，滤波器14'),xlabel('频率'),ylabel('功率谱db')

（5）.最大似然法： 思想：最大似然法原理是让信号通过一个滤波器，选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过，同时，使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。在这个条件下，信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。可以证明,如果信号x 是由一个确定性信号S 加上一个高斯白噪声n 所组成，则上述滤波器的输出是信号S 的最大似然估值。如果n 不是高斯噪声，则上述滤波器的输出是信号S 的最小方差的线性的无偏估值，而且它能得到比使用固定的窗口函数的周期图法更高的分辨率。

程序： fs=20000n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(4) subplot(3,1,1) plot(n,xn)axis([0 0.01 -5 5]) %估计自相关函数 m=-500:500[r,lag]=xcorr(xn,500,'biased')R=[r(501) r(502) r(503) r(504)r(500) r(501) r(502) r(503)r(499) r(500) r(501) r(502)r(498) r(499) r(500) r(501)][V,D]=eig(R)V3=[V(1,3),V(2,3),V(3,3),V(4,3)].'V3=[V(1,4),V(2,4),V(3,4),V(4,4)].'p=0:3wm=[0:0.002*pi:2*pi]B=[(exp(-li)).^(wm'*p)]%li 就是虚数单位j A=B%最小方差功率谱估计 z=A*inv(R)*A'Z=diag(z')pmv=1./Zsubplot(2,1,2) plot(wm/pi,pmv)title('基于最大似然的功率谱估计') ylabel('功率谱幅度（db）') xlabel('角度频率w/pi')

自己分下段就能运行了。

---