矩阵中正弦波的总功率

Question

我正在尝试使用pwelch或periodogram来获取不同频率的信号功率。作为测试，我尝试了一个正弦：

sqrt(2)*sin(x)

它应该有一个一的幂。但是，我的脚本在绘图和总功率中显示了不同的值。我做错了什么，为什么光谱取决于N的样本数量？

clc;
clear all;
fs=1e6;
fn=fs/2;
N=2^12;
df=fs/N;
t=0:1/fs:(N-1)/fs;
nu0=250e3;
nu0=df*floor(nu0/df);
a=sqrt(2);
x=a*sin(2*pi*nu0*t);
[px,f]=pwelch(x,[],0,[],fs);
%[px,f]=periodogram(x,[],[],fs);
semilogy(f,px);
disp(['Total power: ',num2str(sum(px))]);

编辑：澄清我的问题： 我使用以下方法获得了预期的结果：

f=0:df:fn-df;    
X=fft(x);
amp=(abs(X)./N).^2;
amp=2*amp(1:N/2);
amp(1)=amp(1)/2;
plot(f,pow2db(amp));

与 pwelch 和周期图有什么区别？

谢谢！

Answer 1

pwelch和periodogram都给出了功率谱密度。要计算功率，您需要对功率谱密度相对于频率进行积分。将积分近似于一个总和，这意味着您需要将sum(px)乘以频率步长：

sum(px) * (f(2)-f(1))

或者使用带有两个输入的trapz以数值方式计算积分：

trapz(f, px)

在您的示例中，pwelch这给出了

>> [px,f]=pwelch(x,[],0,[],fs);
>> sum(px) * (f(2)-f(1))
ans =
1.000000000000001
>> trapz(f, px)
ans =
0.999999999999995

并带有periodogram：

>> [px,f]=periodogram(x,[],[],fs);
>> sum(px) * (f(2)-f(1))
ans =
1.000000000000000
>> trapz(f, px)
ans =
1.000000000000000