아래 그림과 같이 내가 사인 기능을 다음에 신호를 생성하는, 난 작은 프로그램을 작성 테스트하려면 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를 계산하기 위해 노력하고 있지만 해요 얻은 결과 :
어느 한 사전에 위의 접근 방식
감사 오류를 추적에서 도와 줄 수 하르트 무트 응답 I'vve 코드를 편집하지만 여전히 같은 결과를 얻었다 * UPDATE 후 :
,584,651 :, 입력 된 데이터는 다음과 같다
UPDATE 샘플 여기에 2048의 창 크기를 frequencyand 증가하는 것은 내가 알아 낸 후 : UPDATE 를 결과 창을 사용하여 같은 모양을 여기에 추가 기능을 사용 후 :
첫 번째 프로그램에서는 '시간'을 선언하고 'zeit'을 사용하는 것처럼 보입니다. 컴파일러가 독일어를 아는 경우를 제외하고는 컴파일하지 말아야합니다. – Jakub
아니요 그냥 SO 변수 이름을 바꿉니다 – Engine
@Engine "sinvalues"를 읽기 모드 (rb)에서 열고 fwrite를하고 있습니다. fopen을하는 동안 "wb"를 사용하면 안됩니까? – Icarus3