Tracciare lo spettro di frequenza con C++

Question

Si prega di vedere le modifiche nella risposta sotto questa domanda.

Ho scritto uno script per tracciare lo spettro di frequenza di un segnale sinusoidale con C++. Ecco i passaggi

1. finestra Hanning
2. Applicazione Applicare FFT utilizzando fftw3 biblioteca

Ho tre grafici: Signal, segnale quando viene moltiplicato per la funzione Hanning, e lo spettro di frequenza. Lo spettro di frequenza sembra sbagliato. Dovrebbe avere un picco a 50 Hz. Qualsiasi suggerimento sarebbe apprezzato. Ecco il codice:

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <fftw3.h> 
#include <iostream> 
#include <cmath> 
#include <fstream> 
using namespace std; 

int main() 
{ 
int i; 
double y; 
int N=50; 
double Fs=1000;//sampling frequency 
double T=1/Fs;//sample time 
double f=50;//frequency 
double *in; 
fftw_complex *out; 
double t[N];//time vector 
double ff[N]; 
fftw_plan plan_forward; 

in = (double*) fftw_malloc(sizeof(double) * N); 
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 

for (int i=0; i< N;i++) 
{ 
    t[i]=i*T; 
    ff[i]=1/t[i]; 
    in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform 
    double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window 
    in[i] = multiplier * in[i]; 
    } 

    plan_forward = fftw_plan_dft_r2c_1d (N, in, out, FFTW_ESTIMATE); 

    fftw_execute (plan_forward); 

    double v[N]; 

    for (int i = 0; i < N; i++) 
    { 

    v[i]=20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])/N/2);//Here I have calculated the y axis of the spectrum in dB 

    } 

    fstream myfile; 

    myfile.open("example2.txt",fstream::out); 

    myfile << "plot '-' using 1:2" << std::endl; 

    for(i = 0; i < N; ++i) 

    { 

     myfile << ff[i]<< " " << v[i]<< std::endl; 

    } 

myfile.close(); 

fftw_destroy_plan (plan_forward); 
fftw_free (in); 
fftw_free (out); 
return 0; 
    } 

devo aggiungere che ho tracciato i grafici con gnuplot dopo aver inserito i risultati in example2.txt. Quindi ff [i] vs v [i] dovrebbe darmi lo spettro di frequenze.

Qui ci sono delle piazzole: Spectrum Frequenza e finestra temporale sinusoidale rispettivamente:

Answer 1

Gli intervalli di frequenza erano completamente errati. Secondo http://www.ni.com/white-paper/3995/en/#toc1; la gamma di frequenza e la risoluzione su x -assenti dipendono dalla frequenza di campionamento e N. L'ultimo punto su l'asse delle frequenze dovrebbe essere Fs/2-Fs/N e la risoluzione DF = FS/N .Così ho cambiato il mio script per: (dato che la risoluzione di frequenza è Fs/N come si aumenta il numero di smaples N (o diminuire la frequenza di campionamento Fs) si ottiene risoluzione di frequenza più piccolo e migliori risultati.)

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <fftw3.h> 
#include <iostream> 
#include <cmath> 
#include <fstream> 
using namespace std; 

int main() 
{ 
int i; 
double y; 
int N=550;//Number of points acquired inside the window 
double Fs=200;//sampling frequency 
double dF=Fs/N; 
double T=1/Fs;//sample time 
double f=50;//frequency 
double *in; 
fftw_complex *out; 
double t[N];//time vector 
double ff[N]; 
fftw_plan plan_forward; 

in = (double*) fftw_malloc(sizeof(double) * N); 
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 

for (int i=0; i<= N;i++) 
{ 
t[i]=i*T; 

in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform 
double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window 
in[i] = multiplier * in[i]; 
} 

for (int i=0; i<= ((N/2)-1);i++) 
{ff[i]=Fs*i/N; 
} 
plan_forward = fftw_plan_dft_r2c_1d (N, in, out, FFTW_ESTIMATE); 

fftw_execute (plan_forward); 

double v[N]; 

for (int i = 0; i<= ((N/2)-1); i++) 
{ 

v[i]=(20*log(sqrt(out[i][0]*out[i][0]+ out[i][1]*out[i][1])))/N; //Here I have calculated the y axis of the spectrum in dB 

    } 

fstream myfile; 

myfile.open("example2.txt",fstream::out); 

myfile << "plot '-' using 1:2" << std::endl; 

for(i = 0;i< ((N/2)-1); i++) 

{ 

myfile << ff[i]<< " " << v[i]<< std::endl; 

} 

myfile.close(); 

fftw_destroy_plan (plan_forward); 
fftw_free (in); 
fftw_free (out); 
return 0; 
} 

Answer 2

penso che non si può avere abbastanza campioni, in particolare, riferimento a questo post Electronics.StackExhcange: https://electronics.stackexchange.com/q/12407/84272.

Stai campionando per 50 campioni, quindi 25 scomparti FFT. Stai campionando a 1000 Hz, quindi 1000/2/25 == 250 Hz per i raccoglitori FFT. La risoluzione del tuo cestino è troppo bassa.

Penso che sia necessario ridurre la frequenza di campionamento o aumentare il numero di campioni.

Answer 3

Dal momento che la vostra domanda in su SO, il codice potrebbe usare un po 'di indentazione e stile miglioramento renderlo più facile da leggere.

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <fftw3.h> 
#include <iostream> 
#include <cmath> 
#include <fstream> 
using namespace std; 

int main(){ 
    // use meaningful names for all the variables 
    int i; 
    double y; 
    int N = 550; // number of points acquired inside the window 
    double Fs = 200; // sampling frequency 
    double dF = Fs/N; 
    double T = 1/Fs; // sample time 
    double f = 50; // frequency 
    double *in; 
    fftw_complex *out; 
    double t[N]; // time vector 
    double ff[N]; 
    fftw_plan plan_forward; 

    in = (double*) fftw_malloc(sizeof(double) * N); 
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 

    for (int i = 0; i <= N; i++){ 
     t[i]=i*T; 
     in[i] = 0.7 * sin(2 * M_PI * f * t[i]); // generate sine waveform 
     double multiplier = 0.5 * (1 - cos(2 * M_PI * i/(N-1))); // Hanning Window 
     in[i] = multiplier * in[i]; 
    } 

    for(int i = 0; i <= ((N/2)-1); i++){ 
     ff[i] = (Fs * i)/N; 
    } 

    plan_forward = fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE); 

    fftw_execute(plan_forward); 

    double v[N]; 
    // Here I have calculated the y axis of the spectrum in dB 
    for(int i = 0; i <= ((N/2)-1); i++){ 
     v[i] = (20 * log(sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1])))/N; 
    } 

    fstream myfile; 
    myfile.open("example2.txt", fstream::out); 
    myfile << "plot '-' using 1:2" << std::endl; 

    for(i = 0; i < ((N/2)-1); i++){ 
     myfile << ff[i] << " " << v[i] << std::endl; 
    } 
    myfile.close(); 

    fftw_destroy_plan(plan_forward); 
    fftw_free(in); 
    fftw_free(out); 

    return 0; 
    } 

Il codice può utilizzare altri commenti, soprattutto prima di loop o chiamate di funzione per specificare il valore di ingresso (scopo) e/o il valore di ritorno (risultato).