回复 #14 yangzj 的帖子
那他后面的加窗、STFT的公式怎么回事？而且一个语音信号，就如“开门”两个字的信号也不止10~30ms长度，那该提取哪一段用于做功率谱呢？
回复 #16 zhlong 的帖子
仔细看了下，确实有这个问题，资料中对这个功率谱的得来说的不清楚。这样看起来的话是不是做短时分析后再取平均值。
回复 #17 yangzj 的帖子
^_^，看样子只有写那个资料的人最清楚了。平均一下不知道表达的是什么意思？和直接FFT的区别呢？ 我也是怎么都不懂资料中的方法
这是我原来的程序
function = sheji(file)
data=wavread(file);%语音文件读取
dat(1,:)=data(:,1);%取单路数据
ave=mean(dat);%计算直流分量
dat=dat-ave;%减去直流分量
w=hamming(256);%hamming 窗函数
dat=conv(dat,w);%卷积滤波
n=length(dat);%求dat长度
X=fft(dat);%进行FFT变换
%subplot(2,1,1);
%plot(abs(X));%画图
Y(:,floor(1:n/5))=X(:,floor(1:n/5));%取能量较集中的低频段
%subplot(2,1,2);
Y=abs(Y);
Y=Y.^2;%Y矩阵元素平方
plot(Y);%画图
这个明显我把信号全部都当成了平稳的做
而且我得出的东西就是尖尖的图 我怀疑这个图没有任何意义 有没有可能是用PSD画的 要怎么画呢 从这个资料附录的程序里找到他用的是
=spectrum(x,512,128,hamming(216),22050);
spectrum是老版本的用法，我还没找到它的用法，呵呵 你从哪里看到的？spectrum不是你那种用法啊？
help spectrum
SPECTRUM Spectral Estimation.
H = SPECTRUM.<ESTIMATOR> returns a spectral estimator in H of type
specified by ESTIMATOR.
There are three types of estimators: Power Spectral Density (PSD),
Mean-square Spectrum (MSS), and Pseudo Spectrum.Valid values for
<ESTIMATOR> for each of these three types are listed below. Type
"help spectrum.<ESTIMATOR>" to get help on a specific spectrum
estimator -e.g., "help spectrum.welch"):
----------------------------
Power Spectral Density (PSD)
----------------------------
burg - Burg periodogram - Periodogram
cov- Covariance welch - Welch
mcov - Modified Covariance yulear - Yule-Walker AR
mtm- Thomson multitaper method (MTM)
To calculate the PSD you must create one of the estimators above and
pass it to one of the functions below along with the data.
Fs = 1000; t = 0:1/Fs:.296;
x = cos(2*pi*t*200)+randn(size(t));
h = spectrum.welch; % Create a Welch spectral estimator.
Hpsd = psd(h,x,'Fs',Fs); % Calculate the PSD
plot(Hpsd) % Plot the PSD. 我在网上搜了楼主给的资料，里面有附程序，那里spectrum用的是老版用法
function y=psdfun(x)%函数用来求信号的特征频率向量
=spectrum(x,512,128,hamming(216),22050);
m=0;
for k=2:length(p)-1;
if p(k-1)<p(k)&p(k)>p(k+1)
m=m+1;
y(m)=f(k);
end
end
function frez=eigfrez(x,N)%函数从某人的一组语音信号的特征频率向量中提取出标准特征向量frez，x为这组特征频率向量的细胞矩阵，N为信号数。
bin=zeros(N,60);n=zeros(N,60);y=zeros(N,60);
for i=1:N;
for j=1:length(x{i});
temp=fix(x{i}(1,j)/100);
switch temp
case 0
bin(i,1)=1;n(i,1)=n(i,1)+1;y(i,1)=[(n(i,1)-1)*y(i,1)+x{i}(1,j)]/n(i,1);
case 1
bin(i,2)=1;n(i,2)=n(i,2)+1;y(i,2)=[(n(i,2)-1)*y(i,2)+x{i}(1,j)]/n(i,2);
case 2
bin(i,3)=1;n(i,3)=n(i,3)+1;y(i,3)=[(n(i,3)-1)*y(i,3)+x{i}(1,j)]/n(i,3);
case 3
bin(i,4)=1;n(i,4)=n(i,4)+1;y(i,4)=[(n(i,4)-1)*y(i,4)+x{i}(1,j)]/n(i,4);
case 4
bin(i,5)=1;n(i,5)=n(i,5)+1;y(i,5)=[(n(i,5)-1)*y(i,5)+x{i}(1,j)]/n(i,5);
case 5
bin(i,6)=1;n(i,6)=n(i,6)+1;y(i,6)=[(n(i,6)-1)*y(i,6)+x{i}(1,j)]/n(i,6);
case 6
bin(i,7)=1;n(i,7)=n(i,7)+1;y(i,7)=[(n(i,7)-1)*y(i,7)+x{i}(1,j)]/n(i,7);
case 7
bin(i,8)=1;n(i,8)=n(i,8)+1;y(i,8)=[(n(i,8)-1)*y(i,8)+x{i}(1,j)]/n(i,8);
case 8
bin(i,9)=1;n(i,9)=n(i,9)+1;y(i,9)=[(n(i,9)-1)*y(i,9)+x{i}(1,j)]/n(i,9);
case 9
bin(i,10)=1;n(i,10)=n(i,10)+1;y(i,10)=[(n(i,10)-1)*y(i,10)+x{i}(1,j)]/n(i,10);
case 10
bin(i,11)=1;n(i,11)=n(i,11)+1;y(i,11)=[(n(i,11)-1)*y(i,11)+x{i}(1,j)]/n(i,11);
case 11
bin(i,12)=1;n(i,12)=n(i,12)+1;y(i,12)=[(n(i,12)-1)*y(i,12)+x{i}(1,j)]/n(i,12);
case 12
bin(i,13)=1;n(i,13)=n(i,13)+1;y(i,13)=[(n(i,13)-1)*y(i,13)+x{i}(1,j)]/n(i,13);
case 13
bin(i,14)=1;n(i,14)=n(i,14)+1;y(i,14)=[(n(i,14)-1)*y(i,14)+x{i}(1,j)]/n(i,14);
case 14
bin(i,15)=1;n(i,15)=n(i,15)+1;y(i,15)=[(n(i,15)-1)*y(i,15)+x{i}(1,j)]/n(i,15);
case 15
bin(i,16)=1;n(i,16)=n(i,16)+1;y(i,16)=[(n(i,16)-1)*y(i,16)+x{i}(1,j)]/n(i,16);
case 16
bin(i,17)=1;n(i,17)=n(i,17)+1;y(i,17)=[(n(i,17)-1)*y(i,17)+x{i}(1,j)]/n(i,17);
case 17
bin(i,18)=1;n(i,18)=n(i,18)+1;y(i,18)=[(n(i,18)-1)*y(i,18)+x{i}(1,j)]/n(i,18);
case 18
bin(i,19)=1;n(i,19)=n(i,19)+1;y(i,19)=[(n(i,19)-1)*y(i,19)+x{i}(1,j)]/n(i,19);
case 19
bin(i,20)=1;n(i,20)=n(i,20)+1;y(i,20)=[(n(i,20)-1)*y(i,20)+x{i}(1,j)]/n(i,20);
case 20
bin(i,21)=1;n(i,21)=n(i,21)+1;y(i,21)=[(n(i,21)-1)*y(i,21)+x{i}(1,j)]/n(i,21);
case 21
bin(i,22)=1;n(i,22)=n(i,22)+1;y(i,22)=[(n(i,22)-1)*y(i,22)+x{i}(1,j)]/n(i,22);
case 22
bin(i,23)=1;n(i,23)=n(i,23)+1;y(i,23)=[(n(i,23)-1)*y(i,23)+x{i}(1,j)]/n(i,23);
case 23
bin(i,24)=1;n(i,24)=n(i,24)+1;y(i,24)=[(n(i,24)-1)*y(i,24)+x{i}(1,j)]/n(i,24);
case 24
bin(i,25)=1;n(i,25)=n(i,25)+1;y(i,25)=[(n(i,25)-1)*y(i,25)+x{i}(1,j)]/n(i,25);
case 25
bin(i,26)=1;n(i,26)=n(i,26)+1;y(i,26)=[(n(i,26)-1)*y(i,26)+x{i}(1,j)]/n(i,26);
case 26
bin(i,27)=1;n(i,27)=n(i,27)+1;y(i,27)=[(n(i,27)-1)*y(i,27)+x{i}(1,j)]/n(i,27);
case 27
bin(i,28)=1;n(i,28)=n(i,28)+1;y(i,28)=[(n(i,28)-1)*y(i,28)+x{i}(1,j)]/n(i,28);
case 28
bin(i,29)=1;n(i,29)=n(i,29)+1;y(i,29)=[(n(i,29)-1)*y(i,29)+x{i}(1,j)]/n(i,29);
case 29
bin(i,30)=1;n(i,30)=n(i,30)+1;y(i,30)=[(n(i,30)-1)*y(i,30)+x{i}(1,j)]/n(i,30);
case 30
bin(i,31)=1;n(i,31)=n(i,31)+1;y(i,31)=[(n(i,31)-1)*y(i,31)+x{i}(1,j)]/n(i,31);
case 31
bin(i,32)=1;n(i,32)=n(i,32)+1;y(i,32)=[(n(i,32)-1)*y(i,32)+x{i}(1,j)]/n(i,32);
case 32
bin(i,33)=1;n(i,33)=n(i,33)+1;y(i,33)=[(n(i,33)-1)*y(i,33)+x{i}(1,j)]/n(i,33);
case 33
bin(i,34)=1;n(i,34)=n(i,34)+1;y(i,34)=[(n(i,34)-1)*y(i,34)+x{i}(1,j)]/n(i,34);
case 34
bin(i,35)=1;n(i,35)=n(i,35)+1;y(i,35)=[(n(i,35)-1)*y(i,35)+x{i}(1,j)]/n(i,35);
case 35
bin(i,36)=1;n(i,36)=n(i,36)+1;y(i,36)=[(n(i,36)-1)*y(i,36)+x{i}(1,j)]/n(i,36);
case 36
bin(i,37)=1;n(i,37)=n(i,37)+1;y(i,37)=[(n(i,37)-1)*y(i,37)+x{i}(1,j)]/n(i,37);
case 37
bin(i,38)=1;n(i,38)=n(i,38)+1;y(i,38)=[(n(i,38)-1)*y(i,38)+x{i}(1,j)]/n(i,38);
case 38
bin(i,39)=1;n(i,39)=n(i,39)+1;y(i,39)=[(n(i,39)-1)*y(i,39)+x{i}(1,j)]/n(i,39);
case 39
bin(i,40)=1;n(i,40)=n(i,40)+1;y(i,40)=[(n(i,40)-1)*y(i,40)+x{i}(1,j)]/n(i,40);
case40
bin(i,41)=1;n(i,41)=n(i,41)+1;y(i,41)=[(n(i,41)-1)*y(i,41)+x{i}(1,j)]/n(i,41);
case 41
bin(i,42)=1;n(i,42)=n(i,42)+1;y(i,42)=[(n(i,42)-1)*y(i,42)+x{i}(1,j)]/n(i,42);
case 42
bin(i,43)=1;n(i,43)=n(i,43)+1;y(i,43)=[(n(i,43)-1)*y(i,43)+x{i}(1,j)]/n(i,43);
case 43
bin(i,44)=1;n(i,44)=n(i,44)+1;y(i,44)=[(n(i,44)-1)*y(i,44)+x{i}(1,j)]/n(i,44);
case 44
bin(i,45)=1;n(i,45)=n(i,45)+1;y(i,45)=[(n(i,45)-1)*y(i,45)+x{i}(1,j)]/n(i,45);
case 45
bin(i,46)=1;n(i,46)=n(i,46)+1;y(i,46)=[(n(i,46)-1)*y(i,46)+x{i}(1,j)]/n(i,46);
case 46
bin(i,47)=1;n(i,47)=n(i,47)+1;y(i,47)=[(n(i,47)-1)*y(i,47)+x{i}(1,j)]/n(i,47);
case 47
bin(i,48)=1;n(i,48)=n(i,48)+1;y(i,48)=[(n(i,48)-1)*y(i,48)+x{i}(1,j)]/n(i,48);
case 48
bin(i,49)=1;n(i,49)=n(i,49)+1;y(i,49)=[(n(i,49)-1)*y(i,49)+x{i}(1,j)]/n(i,49);
case 49
bin(i,50)=1;n(i,50)=n(i,50)+1;y(i,50)=[(n(i,50)-1)*y(i,50)+x{i}(1,j)]/n(i,50);
case 50
bin(i,51)=1;n(i,51)=n(i,51)+1;y(i,51)=[(n(i,51)-1)*y(i,51)+x{i}(1,j)]/n(i,51);
case 51
bin(i,52)=1;n(i,52)=n(i,52)+1;y(i,52)=[(n(i,52)-1)*y(i,52)+x{i}(1,j)]/n(i,52);
case 52
bin(i,53)=1;n(i,53)=n(i,53)+1;y(i,53)=[(n(i,53)-1)*y(i,53)+x{i}(1,j)]/n(i,53);
case 53
bin(i,54)=1;n(i,54)=n(i,54)+1;y(i,54)=[(n(i,54)-1)*y(i,54)+x{i}(1,j)]/n(i,54);
case 54
bin(i,55)=1;n(i,55)=n(i,55)+1;y(i,55)=[(n(i,55)-1)*y(i,55)+x{i}(1,j)]/n(i,55);
case 55
bin(i,56)=1;n(i,56)=n(i,56)+1;y(i,56)=[(n(i,56)-1)*y(i,56)+x{i}(1,j)]/n(i,56);
case 56
bin(i,57)=1;n(i,57)=n(i,57)+1;y(i,57)=[(n(i,57)-1)*y(i,57)+x{i}(1,j)]/n(i,57);
case 57
bin(i,58)=1;n(i,58)=n(i,58)+1;y(i,58)=[(n(i,58)-1)*y(i,58)+x{i}(1,j)]/n(i,58);
case 58
bin(i,59)=1;n(i,59)=n(i,59)+1;y(i,59)=[(n(i,59)-1)*y(i,59)+x{i}(1,j)]/n(i,59);
case 59
bin(i,60)=1;n(i,60)=n(i,60)+1;y(i,60)=[(n(i,60)-1)*y(i,60)+x{i}(1,j)]/n(i,60);
end
end
end
t=sum(bin);m=0;
for k=1:60;
if t(1,k)>0;
m=m+1;
frez(m)=sum(y(:,k))/t(1,k);
end
end
function d=distance(from,to)%函数用来求特征向量from与标准向量to的距离d
Nf=length(from);
Nt=length(to);
for i=1:Nf;
iffrom(i)<6000;
j=j+1;
dis(j)=min((to-from(i)*ones(1,Nt)).^2);
end
end
d=sum(dis)/(j*10^4); 辛苦yangzj版主了 没找到这个函数的老版用法的说明。
调试运行了下，=spectrum(x,nfft,noverlap,window,fs);
各参数的含义：x为分析序列；
nfft为FFT点数；
noverlap是重叠采样点数；
window为窗，长度可自定
返回功率谱
做法还是重叠分段，分别求功率谱，再做平均
[ 本帖最后由 yangzj 于 2007-5-26 21:49 编辑 ] 是不是直接运行上面的就可以了？
我在台湾网上找到了这个函数FFTONESIDE
这是他的使用例子
% 顯示一個語音音框的單邊頻譜
=wavread('welcome.wav');
signal=y(2047:2047+237-1);
=fftOneSide(signal, fs, 1);
斑竹大人看看啊 这个得出的中间第二个图可以不？
M函数见附件
[ 本帖最后由 lovelydeath 于 2007-5-26 22:16 编辑 ]
回复 #26 lovelydeath 的帖子
这个只是做一帧的单边谱，如果按你上面文章里程序的做法的话，他是做多帧（帧间可有重叠）再做平均 你可参照我上面对spectrum的说明来做：令index=1:NWindow;设信号为x;
1、令y＝x(index),并对y加长度为NWindow的窗；
2、做yf=fft(y,NFFT);
3、得到y的单边功率谱: yP=(abs(yf(1:NFFT/2))).^2/(NFFT);
4、令index=index+NWindow－Noverlap（Noverlap为重叠采样点数）
重复1－4，直到到x终点；
最后把所有的yP做平均。
大概是这样，具体你实现下