//////////////////////////////////////////////////////////
// 白噪声中单频信号功率谱估计 //
//////////////////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#define PI 3.14159265
#define FRAME 512
#define VARIANCE 0.1 //AWGN方差
double Freq[3] = {0.05,0.15,0.3};
double Amp[3] = {2.0, 1.0, 3.0};
double SinData[FRAME];
double WhiteNoise[FRAME];
double InputData[FRAME];
double HammingWin[FRAME];
double PwrSpect[FRAME],LogSpect[FRAME];
//产生正弦波信号
void F_CreatSin(void)
{
int i;
double LPhase[3];
//产生随机相位
srand((unsigned)time(NULL));
LPhase[0] = 2.0 * PI * ((double)rand()/(double)RAND_MAX);
LPhase[1] = 2.0 * PI * ((double)rand()/(double)RAND_MAX);
LPhase[2] = 2.0 * PI * ((double)rand()/(double)RAND_MAX);
for(i=0;i<FRAME;i++)
{
SinData[i] = Amp[0] * sin(2.0 * PI * Freq[0] * ((double)i) + LPhase[0])
+ Amp[1] * sin(2.0 * PI * Freq[1] * ((double)i) + LPhase[1])
+ Amp[2] * sin(2.0 * PI * Freq[2] * ((double)i) + LPhase[2]);
}
}
//白噪声子程序
void F_CreatWhiteNoise(double fSigma)
{
int i;
int nTemp;
double RandomData1[FRAME],RandomData2[FRAME];
srand((unsigned)time(NULL));
for(i=0;i<FRAME;i++)
{
do
{
nTemp = rand();
}
while((nTemp==0)||(nTemp==RAND_MAX));
RandomData1[i]=(double)nTemp/(double)RAND_MAX;
}
for(i=0;i<FRAME;i++)
{
do
{
nTemp = rand();
}
while((nTemp==0)||(nTemp==RAND_MAX));
RandomData2[i]=(double)nTemp/(double)RAND_MAX;
}
for(i=0;i<FRAME;i++)
{
WhiteNoise[i] = sqrt(-2.0*fSigma*log(RandomData1[i]))*cos(2.0*PI*RandomData2[i]);
}
}
//正弦信号中加噪声
void F_AddNoise(void)
{
int i;
for(i=0;i<FRAME;i++)
{
InputData[i] = SinData[i] + WhiteNoise[i];
}
}
//利用DFT计算功率谱
void F_EstimateSpectum(void)
{
int i,nLoop;
double LWinDada[FRAME];
double LSumRe,LSumIm;
//生成Hamming窗函数
for(i=0;i<FRAME;i++)
{
HammingWin[i] = 0.54 + 0.46 * cos((double)(i-(FRAME/2))*2.0*PI/(double)FRAME);
}
//数据加窗处理
for(i=0;i<FRAME;i++)
{
LWinDada[i] = InputData[i] * HammingWin[i];
}
//功率谱估计
for(nLoop=0;nLoop<FRAME;nLoop++)
{
LSumRe = LSumIm = 0.0;
for(i=0;i<FRAME;i++)
{
LSumRe += LWinDada[i] * cos(2.0*PI*double(i*nLoop)/(double)FRAME);
LSumIm += LWinDada[i] * sin(2.0*PI*double(i*nLoop)/(double)FRAME);
}
PwrSpect[nLoop] = (LSumRe*LSumRe + LSumIm*LSumIm)/double(FRAME);
LogSpect[nLoop] = 10.0*log10(PwrSpect[nLoop]);
}
}
////////////////////////////////////////////////////////////////////
void main(void)
{
FILE *fp1,*fp2,*fp3;
int i;
F_CreatSin();
F_CreatWhiteNoise(VARIANCE);
F_AddNoise();
F_EstimateSpectum();
//存储输入数据
if((fp1 = fopen("input.dat","w+")) == NULL)
{
printf("Open file input.dat error !");
exit(1);
};
for(i=0;i<FRAME;i++)
{
fprintf(fp1,"%d %f\n",i,InputData[i]);
}
fclose(fp1);
//存储功率谱参数
if((fp2 = fopen("spectdata.dat","w+")) == NULL)
{
printf("Open file spectdata.dat error !");
exit(1);
};
for(i=0;i<FRAME;i++)
{
fprintf(fp2,"%d %f\n",i,PwrSpect[i]);
}
fclose(fp2);
//存储功率谱参数(对数方式)
if((fp3 = fopen("logspect.dat","w+")) == NULL)
{
printf("Open file logspect.dat error !");
exit(1);
};
for(i=0;i<FRAME;i++)
{
fprintf(fp3,"%d %f\n",i,LogSpect[i]);
}
fclose(fp3);
printf("End\n");
}