PSD使用FFTW Halfcomplex变换
我问了一个类似的问题,这个问题得到了解答,但当我尝试按照自己的方式进行操作时,我会得到“奇怪的”值。 我想得到一个sin波的PSD使用半复数变换,如:
#include #include #include #include #include #define PI 3.141592653589793 int main (){ double* inputValues; double* outputValues; double realVal; double imagVal; double powVal=0.0; double absVal; double timer; fftw_plan plan; double timeIntervall= 1.0; // 1sec int numberOfSamples =512; double timeSteps = timeIntervall/numberOfSamples; float frequency=10.0; float dcValue = 0.2; float value=0.0; int index=0; // allocating the memory for the fftw arrays inputValues = (double*) fftw_malloc(sizeof(double)* numberOfSamples); outputValues = (double *) fftw_malloc(sizeof(double)*(numberOfSamples/*2*/)); plan = fftw_plan_r2r_1d(numberOfSamples,inputValues,outputValues,FFTW_R2HC,FFTW_ESTIMATE); for (timer = 0; timer<=timeIntervall; timer += timeSteps){ value = sin(2*M_PI*frequency*timer) +dcValue; inputValues[index++] = value; } fftw_execute(plan); for (index=0;index<=numberOfSamples/*2*/;index++){ powVal = outputValues[index]*outputValues[index]+outputValues[numberOfSamples-index]*outputValues[numberOfSamples- index]; if(index==0) powVal/=2; powVal/=numberOfSamples; fprintf(stdout," index %d \t PSD value %lf \n",index,powVal); } return 0; }
我得到的价值是:
index 0 PSD value 12.24 // expecting 0.2 ................ ..................... index 10 PSD value 129.99999 // expecting 0.5 ........ ....... index 502 PSD value 127.9999 // expecting 0.5 ...................... ...................... index 512 PSD value 12.24 // expecting 0.2
否则PSD值为零,峰值的位置是正确的,但它们的值不知道为什么?
提前致谢 !
UPDATE
我解决了,但我不明白为什么它有效,所以我不会把它作为答案:这是我在代码中改变了:
....................................... fftw_execute_r2r(plan_r2hc, in, out); powVal = outputValues[0]*outputValues[0]; powVal /= (numberOfSamples*numberOfSamples)/2; ///WHY ?????? index = 0; fprintf(stdout," index %d \t PSD value %lf \t \t %lf \n",index,powVal,outputValues[index]); for (index =1; index<numberOfSamples/2;index++){ powVal = outputValues[index]*outputValues[index]+outputValues[numberOfSamples-index]*outputValues[numberOfSamples- index]; powVal/=(numberOfSamples*numberOfSamples)/2; //WHY????? fprintf(stdout," index %d \t PSD value %lf \t \t %lf \n",index,powVal,outputValues[index]); }
结果是准确的,我希望得到任何解释为什么我应该划分windowsSize和on 2的平方? 再次感谢你的帮助 !
如数值配方中所述,功率谱密度(PSD)的归一化可能会根据您对PSD的确切定义而有所不同。 一种可能的定义使得频谱估计的总和对应于时域函数的均方幅度:
其中P(k)与FFT输出X(k)相关:
对于一些比例因子S
这似乎是您根据预期结果使用的定义。
应用Parseval定理 :
到定义产量:
或S = 1 / N 2 。
这当然假设您使用整个频谱。 另一方面, FFTW的半复数格式 ,顾名思义,只给你一半的频谱(另一半是实值输入对称的)。 通过将该对称谱值的等值平方幅度相加,可以得到频率分量的总功率。 这是你获得2的因子。 请注意,实值索引0和numberOfSamples/2
没有对应的对称输出,因此在这些情况下不应乘以2。
因此,您的PSD应计算为:
powVal = outputValues[0]*outputValues[0]; powVal /= (numberOfSamples*numberOfSamples); for (index =1; index
注意:直流分量的均方幅度应为0.2 * 0.2 = 0.04(不是0.2,表示为指数0的预期值)。