怎样求功率谱？我编的程序怎么不对，大家帮忙看看

heiweibo

计算原始信号的功率谱
把声频信号s(i)分帧，每帧长度为512。采样频率为44.1Hz,则信号功率谱密度的计算公

式为10*log(abs(h(i).*s(i)'.*exp(-j*2*pi*k*n/N)/N)^2);

在做这个算法的时候我的程序是这样的，但是结果是不对，大家帮忙看看吧，是那里的问题，谢谢了
clear  all;
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

songzy41

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 编辑 ]

yuanlive

问一个相关的问题，我在做能量普密度的时候，为什么最后得出的都是负值？这样合理吗？
[audio_data,fs,nbits]=wavread('d:/audio/test48_double.wav');
siz=wavread('d:/audio/test48_double.wav','size');
w=audio_data';
s=w([1:1],:);
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

songzy41

正如yangzj在另一帖子中对你的程序所说，由于取了对数，功率谱密度有可能是负值。同时在你的运算中用DFT，而实际上可用FFT，运算速度快了很多，精度一样。