Plotten Frequenzspektrum mit C++

Question

Bitte beachten Sie die Edits in der Antwort unter dieser Frage.

Ich habe ein Skript geschrieben, um das Frequenzspektrum eines Sinussignals mit C++ zu zeichnen. Hier sind die Schritte

1. Anwendung Hanning Fenster
2. Nehmen FFT fftw3 Bibliothek

Ich habe drei Graphen: Signal, Signal, wenn auf Hanning-Funktion multipliziert wird, und das Frequenzspektrum. Das Frequenzspektrum sieht falsch aus. Es sollte eine Spitze bei 50 Hz haben. Jeder Vorschlag wäre willkommen. Hier ist der Code:

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

Ich muss hinzufügen, dass ich die Graphen aufgetragen haben mit gnuplot nach Erhalt der Ergebnisse in example2.txt einsetzen. Also sollte ff [i] vs v [i] mir das Frequenzspektrum geben.

Hier sind die Plots: Frequenzspektrum und Sinuszeitfenster jeweils:

Answer 1

Meine Frequenzintervalle waren völlig falsch. Nach http://www.ni.com/white-paper/3995/en/#toc1; Der Frequenzbereich und die Auflösung auf der x -Achse hängen von der Abtastrate und N ab. Der letzte Punkt auf der Frequenzachse sollte Fs/2-Fs/N und die Auflösung dF = FS/N sein. Also habe ich mein Skript geändert in: (seit Frequenzauflösung ist Fs/N wie Sie erhöhen die Anzahl der smaples N (oder Abtastfrequenz Fs verringern) Sie kleinere Frequenzauflösung erhalten und bessere Ergebnisse.)

#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

Ich glaube, Sie nicht genug Proben haben können, vor allem, verweist auf diesen Electronics.StackExhcange Beitrag: https://electronics.stackexchange.com/q/12407/84272.

Sie sampling für 50 Proben, also 25 FFT-Bins. Sie nehmen eine Abtastung mit 1000 Hz vor, also 1000/2/25 = 250 Hz pro FFT-Bins. Ihre bin-Auflösung ist zu niedrig.

Ich denke, Sie müssen die Abtastfrequenz senken oder die Anzahl der Proben erhöhen.

Answer 3

Da Ihre Frage in auf SO, könnte Ihr Code eine Vertiefung und Stil Verbesserung verwenden, um mache es einfacher zu lesen.

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

Ihr Code kann mehr Kommentare verwenden, besonders vor dem Schleifen oder Funktionsaufrufe dessen Eingangswert (Ziel) und/oder die Rückkehr Wert (Ergebnis) zu spezifizieren.