怎样求功率谱？我编的程序怎么不对，大家帮忙看看
计算原始信号的功率谱把声频信号s(i)分帧，每帧长度为512。采样频率为44.1Hz,则信号功率谱密度的计算公
式为10*log(abs(h(i).*s(i)'.*exp(-j*2*pi*k*n/N)/N)^2);
在做这个算法的时候我的程序是这样的，但是结果是不对，大家帮忙看看吧，是那里的问题，谢谢了
clearall;
close all;
sp=wavread('a');
N=512;T=fix(length(sp)/N);%分帧
s=zeros(1,N)；
h=hanning(N);
sum=0;j=sqrt(-1);
for p=1:T
s=sp((p-1)*N+1:p*N);
for k=1:N
for i=1:N
sum=sum+h(i)*s(i)*exp(-j*2*pi*k*i/N);
end
sum=abs(sum);
x(k)=20*log10(sum);
end
result(:,p)=x';
end 1，在程序中主要的错是sum=0放错了位置，应把该语句放在二层for之间：
for k=1:N
sum=0;
for i=1:N
而在公式中Σ以后的算式实际上就是计算DFT（h*s的DFT），因此计算速度特慢。
2，程序中的DFT运算可用FFT去替代：
for p=1:T
s=sp((p-1)*N+1:p*N);
x=fft(h.*s);
x=20*log10(abs(x));
result(:,p)=x;
end
这样运算速度可快多了，得到是结果是一样的。
[ 本帖最后由 songzy41 于 2006-9-22 17:54 编辑 ] 问一个相关的问题，我在做能量普密度的时候，为什么最后得出的都是负值？这样合理吗？
=wavread('d:/audio/test48_double.wav');
siz=wavread('d:/audio/test48_double.wav','size');
w=audio_data';
s=w(,:);
N=1024;
h=hanning(N);
j=sqrt(-1);
S=s(1:N);
S=s(1:N);
for k=1:1:512
sum=0;
for m=1:1:N
tmp1=h(m)*S(m)*exp(-j*2*pi*k*m/N);
sum=tmp1+sum;
end
x(k)=(abs(sum/N))^2;
X(k)=10*log10(x(k));
end 正如yangzj在另一帖子中对你的程序所说，由于取了对数，功率谱密度有可能是负值。同时在你的运算中用DFT，而实际上可用FFT，运算速度快了很多，精度一样。
页:
[1]