1
$\begingroup$

I've already asked this question on SO, but didn't get an answer, I tqake a shot here. I am trying to get the PSD of a real data set by making use of fftw3 library To test I wrote a small program as shown below ,that generates the a signal which follows sinusoidal function

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

int main (){
    double  value= 0.0;
    float frequency = 5;
    float  length = 2;
    int index = 4;
    double 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;

}

Now I'm reading the output file of above program and trying to calculate its PSD like as shown below

#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;

}

Iam not sure about the correctness of the way I am doing this, but below are the results i have obtained:

PSD

Can any one help me in tracing the errors of the above approach

and the input data look like :

sinus

Thanks in advance

$\endgroup$
3
  • 1
    $\begingroup$ In your code listing the sampling frequency and the frequency of your tone are not defined very well. zeit is not initialized. Also zeit = n/F, where n is an integer. Then $\sin(2\pi*F*zeit)= \sin(2\pi*F*n/F)=\sin(2\pi*n)=0$, so the output will be a constant. $\endgroup$
    – David
    Jan 7 '15 at 20:38
  • $\begingroup$ your sine production is essentailly broken, as @David says. $\sin(2\pi f \frac nf)\equiv \sin(2\pi n)\equiv 0 \quad \forall n\in\mathbb N$. the fact you're seeing something like a sine at all is only due to the limited machine resolution! $\endgroup$ Jul 31 '16 at 14:56
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because the data analyzed here is produced in a broken way and answering the question hence doesn't solve any relevant problem. $\endgroup$ Jul 31 '16 at 14:57
0
$\begingroup$

Well, first of all you should note that you are taking batches of 512 samples. It basically means that the first FFT will perform a transform of the first 512 samples, slightly more than a quarter of a period, assuming that signal is periodic. (The windowing tries to compensate for this to some extent.) However, the information from that signal is not really what brings you any information about the spectrum.

While I haven't read the code, I suggest that you try one of these things (maybe both):

  • Increase the signal frequency

  • Increase the window length

Ideally, there should be at least one full period of the sinusoidal in the window if you want to interpret the result in a straightforward way.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.