2014-07-11 1 views
1

아래 그림과 같이 내가 사인 기능을 다음에 신호를 생성하는, 난 작은 프로그램을 작성 테스트하려면 fftw3 라이브러리 이용하여 설정된 실제 데이터의 PSD를 얻기 위해 노력하고파워 스펙트럼 밀도를 계산

#include <stdio.h> 
#include <math.h> 
#define PI 3.14 

int main(){ 
    double value= 0.0; 
    float frequency = 5; 
    int i = 0 ; 
    double time = 0.0; 
    FILE* outputFile = NULL; 
    outputFile = fopen("sinvalues","wb+"); 
    if(outputFile==NULL){ 
     printf(" couldn't open the file \n"); 
     return -1; 
    } 

    for (i = 0; i<=5000;i++){ 

     value = sin(2*PI*frequency*zeit); 
     fwrite(&value,sizeof(double),1,outputFile); 
     zeit += (1.0/frequency); 
    } 
    fclose(outputFile); 
    return 0; 

} 

지금은 아래에 있습니다 위의 프로그램의 출력 파일을 읽고

#include <stdio.h> 
#include <fftw3.h> 
#include <complex.h> 
#include <stdlib.h> 
#include <math.h> 
#define PI 3.14 
int main(){ 
    FILE* inp = NULL; 
    FILE* oup = NULL; 
    double* value;// = 0.0; 
    double* result; 
    double spectr = 0.0 ; 
    int windowsSize =512; 
    double power_spectrum = 0.0; 
    fftw_plan plan; 

    int index=0,i ,k; 
    double multiplier =0.0; 
    inp = fopen("1","rb"); 
    oup = fopen("psd","wb+"); 

    value=(double*)malloc(sizeof(double)*windowsSize); 
    result = (double*)malloc(sizeof(double)*(windowsSize)); // what is the length that I have to choose here ? 
     plan =fftw_plan_r2r_1d(windowsSize,value,result,FFTW_R2HC,FFTW_ESTIMATE); 

    while(!feof(inp)){ 

     index =fread(value,sizeof(double),windowsSize,inp); 
      // zero padding 
     if(index != windowsSize){ 
      for(i=index;i<windowsSize;i++){ 
        value[i] = 0.0; 
         } 

     } 


     // windowing Hann 

     for (i=0; i<windowsSize; i++){ 
      multiplier = 0.5*(1-cos(2*PI*i/(windowsSize-1))); 
      value[i] *= multiplier; 
     } 


     fftw_execute(plan); 


     for(i = 0;i<(windowsSize/2 +1) ;i++){ //why only tell the half size of the window 
      power_spectrum = result[i]*result[i] +result[windowsSize/2 +1 -i]*result[windowsSize/2 +1 -i]; 
      printf("%lf \t\t\t %d \n",power_spectrum,i); 
      fprintf(oup," %lf \n ",power_spectrum); 
     } 

    } 
    fclose(oup); 
    fclose(inp); 
    return 0; 

} 

스피 나는이 일을하고있는 방법의 정확성에 대해 확실하지 아래와 같이 같은 자사의 PSD를 계산하기 위해 노력하고 있지만 해요 얻은 결과 :

PSD

어느 한 사전에 위의 접근 방식

감사 오류를 추적에서 도와 줄 수 하르트 무트 응답 I'vve 코드를 편집하지만 여전히 같은 결과를 얻었다 * UPDATE 후 :

,584,651 : PSD2

, 입력 된 데이터는 다음과 같다sinus

UPDATE 샘플 여기에 2048의 창 크기를 frequencyand 증가하는 것은 내가 알아 낸 후 : enter image description here UPDATE 를 결과 창을 사용하여 같은 모양을 여기에 추가 기능을 사용 후 : enter image description here

+0

첫 번째 프로그램에서는 '시간'을 선언하고 'zeit'을 사용하는 것처럼 보입니다. 컴파일러가 독일어를 아는 경우를 제외하고는 컴파일하지 말아야합니다. – Jakub

+0

아니요 그냥 SO 변수 이름을 바꿉니다 – Engine

+0

@Engine "sinvalues"를 읽기 모드 (rb)에서 열고 fwrite를하고 있습니다. fopen을하는 동안 "wb"를 사용하면 안됩니까? – Icarus3

답변

2

잘못된 출력 값을 파워 스펙트럼 라인에 결합합니다.역순으로 끝에 및 windowsSize/2 - 1 허수 값의 시작 부분에 windowsSize/2 + 1 실제 값이 있습니다. 하여 첫 번째 (0Hz에서)의 상상의 구성 요소와 마지막 (나이 퀴 스트 주파수) 스펙트럼 라인은 0

int spectrum_lines = windowsSize/2 + 1; 
power_spectrum = (double *)malloc(sizeof(double) * spectrum_lines); 

power_spectrum[0] = result[0] * result[0]; 
for (i = 1 ; i < windowsSize/2 ; i++) 
    power_spectrum[i] = result[i]*result[i] + result[windowsSize-i]*result[windowsSize-i]; 
power_spectrum[i] = result[i] * result[i]; 

입니다 그리고 작은 실수가 있기 때문이다 : 당신은 입력 신호에 아닌 윈도우 함수를 적용해야 제로 패딩 부로 보낸다.

ADD-ON은 :

테스트 프로그램은 정현파 신호의 5001 개 샘플을 생성하고이 신호의 512 개 샘플을 읽고 분석 할 수 있습니다. 그 결과, 일정 기간 만 분석 할 수 있습니다. 하드 컷오프 신호로 인해 예측할 수없는 에너지 레벨을 가진 넓은 스펙트럼의 에너지를 포함합니다. PI는 사용하지 않고 예측 가능한 계산을 수행하기에는 정확하지 않은 3.41 만 사용하기 때문입니다.

정수의 피리어드가 512 샘플의 분석 창과 정확히 일치해야합니다. 따라서 테스트 신호 작성 프로그램에서 테스트 신호에서 정확하게 numberOfPeriods 마침표를 변경해야합니다 (예 : numberOfPeriods=1은 시노이드의 한주기에 정확히 512 샘플, 2 => 256,3 => 512/3의 마침표가 있음을 의미합니다 , 4 => 128, ...). 이렇게하면 특정 스펙트럼 라인에서 에너지를 생성 할 수 있습니다. windowSize은 크기가 다르기 때문에 이러한 노력이 쓸모 없기 때문에 두 프로그램에서 동일한 값을 가져야합니다.

#define PI 3.141592653589793 // This has to be absolutely exact! 

int windowSize = 512;  // Total number of created samples in the test signal 
int numberOfPeriods = 64; // Total number of sinoid periods in the test signal 
for (n = 0 ; n < windowSize ; ++n) { 
    value = sin((2 * PI * numberOfPeriods * n)/windowSize); 
    fwrite(&value, sizeof(double), 1, outputFile); 
} 
+0

windowsSize는 512 왜 spectral_line이 필요합니까? – Engine

+0

필요 없어. 스펙트럼 선의 수를 절대적으로 명확하게하는 것이 유용합니다. 당신은 그것을 건너 뛸 수 있습니다. 그것은 단지 프로그램의 주석에 "왜 창 크기의 절반 만 말해라"고 물었 기 때문입니다. 그게 이유야. –

+0

코드를 변경 한 후에도 여전히 비슷한 결과를 얻었습니다 – Engine

1

귀하의 질문이 확실하지 않습니다. 제공되는 정보와 함께 귀하의 결과는 합리적으로 보입니다.

PSD는 자기 상관 함수의 푸리에 변환입니다. 사인파 입력의 경우 AC 기능은 주기적이므로 PSD는 사용자가 음모를 꾸미는 것처럼 톤을 갖습니다.

내 '대답'은 디버깅에 대한 기본적인 생각입니다. 방정식을 게시 할 수 있다면 모든 참여자가 더 쉬울 것입니다. 요즘 SE에는 신호 처리 섹션이 있다는 것을 알고있을 것입니다.

먼저 AC 기능의 플롯을 제공해야합니다. 표시 한 PSD의 역 FT는주기적인 톤의 선형 조합입니다.

둘째, 창을 제거하고 상자를 만들거나 가능한 경우 단계를 건너 뜁니다.

셋째, DFT를 FFT로 바꾸어보십시오 (필자는 fftw3 라이브러리 문서 만 다뤘습니다. 아마도 이것은 옵션입니다).

마지막으로 화이트 노이즈를 입력하려고합니다. Bernoulli dist 또는 Gaussian dist를 사용할 수 있습니다. 샘플 AC는 그렇지 않지만 AC는 델타 기능입니다. 이렇게하면 (샘플) 흰색 PSD 분포가 제공됩니다.

이 제안 사항이 도움이되기를 바랍니다.

+0

내가 얻은 결과가 올바르지 않습니다. DC! = 0은 DFT에 관한 것이고 FFTW는 PSF를 계산할 때이 functio를 사용하도록 제안합니다. 결과가 정기적 인 이유는 내가 동등한 봉우리가 아닌 이유입니까? – Engine

+0

AC 플롯을 제공했다면 도움이 될 것입니다. DC 레벨이 0이 아닌 경우이 문은 올바르지 않습니다.백색 잡음 시리즈의 PSD는 주파수 = 0을 포함하여 어디에서나 볼 수 있습니다. – JayInNyc

2

예상 출력 기능에 대한 설명입니다.

  • 귀하의 의견은 순수한 실제 가치가있는 기능입니다. DFT의 결과에는 복잡한 값이 있습니다. 그래서 변수를 double이 아닌 fftw_complex * out으로 선언해야합니다.

  • 일반적으로 입력 값의 수는 출력 값의 수와 같습니다. 그러나 dft의 출력 스펙트럼에는 음의 주파수뿐만 아니라 양수에 대한 복소 진폭이 포함됩니다.

  • 순수한 실수 입력에 대한 특별한 경우에, 양의 주파수의 진폭은 음의 주파수의 진폭의 공액 복소수 값입니다. 그럴 경우 긍정적 인 스펙트럼의 주파수 만 계산됩니다. 이는 복소 출력 값의 수가 실제 입력 값의 수인 의 절반임을 의미합니다.

  • 입력이 단순한 사인파 인 경우 스펙트럼에는 단일 주파수 구성 요소 만 포함됩니다. 이것은 10, 100, 1000 또는 그 이상의 입력 샘플에 해당됩니다. 다른 모든 값은 0입니다. 따라서 엄청난 수의 입력 값으로 작업하는 것은 의미가 없습니다.

  • 입력 데이터 세트에 단일 마침표가 포함 된 경우 복소 출력 값은 out [1]에 포함 된 입니다.입력 데이터 세트가 귀하의 경우 5, M 완료 기간을 포함하면 경우

  • , 는 그래서 결과는 [5]

  • 난 당신의 코드를 일부 수정했다 밖으로에 저장됩니다. 몇 가지 사실을 분명히하기.

    #include <iostream> 
    #include <stdio.h> 
    #include <math.h> 
    #include <complex.h> 
    #include "fftw3.h" 
    
    int performDFT(int nbrOfInputSamples, char *fileName) 
    { 
        int nbrOfOutputSamples; 
        double *in; 
        fftw_complex *out; 
        fftw_plan p; 
    
        // In the case of pure real input data, 
        // the output values of the positive frequencies and the negative frequencies 
        // are conjugated complex values. 
        // This means, that there no need for calculating both. 
        // If you have the complex values for the positive frequencies, 
        // you can calculate the values of the negative frequencies just by 
        // changing the sign of the value's imaginary part 
        // So the number of complex output values (amplitudes of frequency components) 
        // are the half of the number of the real input values (amplitutes in time domain): 
        nbrOfOutputSamples = ceil(nbrOfInputSamples/2.0); 
    
        // Create a plan for a 1D DFT with real input and complex output 
        in = (double*) fftw_malloc(sizeof(double) * nbrOfInputSamples); 
        out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nbrOfOutputSamples); 
        p = fftw_plan_dft_r2c_1d(nbrOfInputSamples, in, out, FFTW_ESTIMATE); 
    
        // Read data from input file to input array 
        FILE* inputFile = NULL; 
        inputFile = fopen(fileName,"r"); 
        if(inputFile==NULL){ 
         fprintf(stdout,"couldn't open the file %s\n", fileName); 
         return -1; 
        } 
        double value; 
        int  idx = 0; 
        while(!feof(inputFile)){ 
         fscanf(inputFile, "%lf", &value); 
         in[idx++] = value;  
        } 
        fclose(inputFile); 
    
        // Perform the dft 
        fftw_execute(p); 
    
        // Print output results 
        char outputFileName[] = "dftvalues.txt"; 
        FILE* outputFile = NULL; 
        outputFile = fopen(outputFileName,"w+"); 
        if(outputFile==NULL){ 
         fprintf(stdout,"couldn't open the file %s\n", outputFileName); 
         return -1; 
        } 
        double realVal; 
        double imagVal; 
        double powVal; 
        double absVal; 
    
        fprintf(stdout, "  Frequency Real  Imag  Abs  Power\n"); 
        for (idx=0; idx<nbrOfOutputSamples; idx++) { 
         realVal = out[idx][0]/nbrOfInputSamples; // Ideed nbrOfInputSamples is correct! 
         imagVal = out[idx][1]/nbrOfInputSamples; // Ideed nbrOfInputSamples is correct! 
         powVal = 2*(realVal*realVal + imagVal*imagVal); 
         absVal = sqrt(powVal/2); 
         if (idx == 0) { 
          powVal /=2; 
         } 
         fprintf(outputFile, "%10i %10.4lf %10.4lf %10.4lf %10.4lf\n", idx, realVal, imagVal, absVal, powVal); 
         fprintf(stdout,  "%10i %10.4lf %10.4lf %10.4lf %10.4lf\n", idx, realVal, imagVal, absVal, powVal); 
         // The total signal power of a frequency is the sum of the power of the posive and the negative frequency line. 
         // Because only the positive spectrum is calculated, the power is multiplied by two. 
         // However, there is only one single line in the prectrum for DC. 
         // This means, the DC value must not be doubled. 
    
        } 
        fclose(outputFile); 
    
    
        // Clean up 
        fftw_destroy_plan(p); 
        fftw_free(in); fftw_free(out); 
    
        return 0; 
    
    } 
    
    int main(int argc, const char * argv[]) { 
        // Set basic parameters 
        float timeIntervall = 1.0;  // in seconds 
        int  nbrOfSamples = 50;  // number of Samples per time intervall, so the unit is S/s 
        double timeStep = timeIntervall/nbrOfSamples; // in seconds 
        float frequency = 5;   // frequency in Hz 
    
        // The period time of the signal is 1/5Hz = 0.2s 
        // The number of samples per period is: nbrOfSamples/frequency = (50S/s)/5Hz = 10S 
        // The number of periods per time intervall is: frequency*timeIntervall = 5Hz*1.0s = (5/s)*1.0s = 5 
    
    
    
        // Open file for writing signal values 
        char fileName[] = "sinvalues.txt"; 
        FILE* outputFile = NULL; 
        outputFile = fopen(fileName,"w+"); 
        if(outputFile==NULL){ 
         fprintf(stdout,"couldn't open the file %s\n", fileName); 
         return -1; 
        } 
    
        // Calculate signal values and write them to file 
        double time; 
        double value; 
        double dcValue = 0.2; 
        int idx = 0; 
        fprintf(stdout, "  SampleNbr Signal value\n"); 
        for (time = 0; time<=timeIntervall; time += timeStep){ 
         value = sin(2*M_PI*frequency*time) + dcValue; 
         fprintf(outputFile, "%lf\n",value); 
         fprintf(stdout, "%10i %15.5f\n",idx++, value); 
        } 
        fclose(outputFile); 
    
        performDFT(nbrOfSamples, fileName); 
        return 0; 
    
    } 
    
    • 의 DFT에 입력 순수한 진짜 인 경우

    출력은 어떤 경우에 복잡하다. 그래서 계획 r2c (RealToComplex)를 사용해야합니다.
  • 신호가 sin (2 * pi * f * t) 인 경우, t = 0에서 시작하여 스펙트럼은 순수한 가상의 f에서 단일 주파수 라인 을 포함합니다.
  • 부호가 sin (2 * pi * f * t + phi)와 같이 위상이 오프셋 인 경우 단일 행의 값은 복소수입니다.
  • 샘플링 주파수가 fs이면 출력 스펙트럼의 범위는 -fs/2 ... + fs/2입니다.
  • 양수 및 음수의 실수 부분은 동일합니다.
  • 양수와 음수의 허수 부에는 반대 부호가 있습니다. 이것을 공액 복합체라고합니다.
  • 양의 스펙트럼의 복잡한 값을 가지고있는 경우 허수 성분의 부호를 변경하여 음수 스펙트럼의 값을 계산할 수 있습니다. 이 때문에 양수 및 음수 sprectrum을 계산할 필요가 없습니다.
  • 한 쪽 밴드는 모든 정보를 보유합니다. 따라서 계획 r2c의 출력 샘플 수는 입력 샘플 수 의 반 + 1입니다.
  • 주파수의 힘을 얻으려면 음의 주파수로 긍정적인 주파수 인 도 고려해야합니다. 그러나 계획 r2c는 스펙트럼의 오른쪽 긍정적 인 절반 만을 제공합니다. 따라서 총 전력을 얻으려면 긍정적 인 측면의 힘을 배가해야합니다.

    fftw3 패키지의 문서는 계획 사용법을 잘 설명합니다. 설명서를 읽으려면 시간을 투자해야합니다.

+0

답장을 보내 주셔서 감사합니다.하지만 몇 가지 질문이 있습니다. 계획이 r2c인데, 입력이 순수한 진짜이므로 r2r를 사용 했으므로 주파수의 절반 만 필요합니다. 2. ouputSamples 크기가 단지 1/2 inputSamples이므로 프로그램이 중단됩니다. 마지막 질문은 왜 powVal에 2를 곱한 것입니까? 도와 주셔서 다시 한 번 감사드립니다! – Engine

+0

@Engine 시스템에서 프로그램이 왜 충돌하는지 알 수 없습니다. OS X Mavericks가있는 Mac을 사용합니다. 이 코드는 Xcode 6.3 베타로 작성되었습니다. 이전 언급에 몇 가지 발언을 추가했습니다. –