Tôi đang cố gắng tạo một chương trình tính PSD của chuỗi thời gian (16384 mẫu), giả sử hình sin ở đây là mã:

// tạo mẫu tội lỗi 

#include 
#include 
#include 

int chính(){
    TẬP TIN *inp =NULL,*inp2;
    giá trị kép = 0,0; 
    tần số thả nổi = 1 // tần số tín hiệu 
    gấp đôi thời gianSec = 1; // thời gian tính bằng giây 
    số int không dấuOFSamples = 2048*8;
    bước kép = timeSec / numberOFSamples; 
    bộ đếm thời gian kép = 0,0;
    float dcValue =0,0;
    chỉ số kép = 0;
     inp = fopen("xoang","wb+");
     inp2=fopen("sinusD","wb+");
    for(timer=0.0 ; bộ đếm thời gian<>
        value= sin(2*M_PI*tần số*bộ đếm thời gian) +dcValue;
        fprintf(inp,"%lf",giá trị);
        fwrite(&value,sizeof(double),1,inp2);

     }

     fclose(inp);
     fclose(inp2);
     return 0;

}

Việc tạo các giá trị sin là chính xác. Bây giờ tôi kiểm tra nó bằng Matlab, PSD phải lớn 1024, đây là mã:

    #include 
#include 
#include 
#include 

#xác định WINDOW_SIZE 1024

int chính(){
  TẬP_TIN* đầu vàoFile = NULL;
  TẬP_TIN* đầu raFile= NULL;
  gấp đôi* dữ liệu đầu vào=NULL; 
  gấp đôi* dữ liệu đầu ra=NULL;
  double* windowData=NULL;
  unsigned int windowSize = 1024;
  int chồng chéo =0;
  int chỉ mục1 =0,index2=0, i=0;
  powVal gấp đôi= 0,0;
  fftw_plan plan_r2hc;
  gấp đôi w[WINDOW_SIZE];
  cho (i=0; i
  w[i] = (1.0 - cos(2.0 * M_PI * i/(WINDOW_SIZE-1))) * 0.5;
  }

// cấp phát bộ nhớ
  inputData = (double*) fftw_malloc(sizeof(double)*windowSize);
  outData= (double*) fftw_malloc(sizeof(double)*windowSize);
  windowData= (double*) fftw_malloc(sizeof(double)*windowSize);
  plan_r2hc = fftw_plan_r2r_1d(windowSize, inputData, windowData, FFTW_R2HC, FFTW_PATIENT);
  //Các tập tin thao tác 
  inputFile = fopen("sinusD","rb");
  outFile= fopen("windowingResult","wb+");
  if(inputFile==NULL){
    printf("Không thể mở được tệp đầu vào hoặc đầu ra \n");
    return -1;
  }

  while((i=fread(inputData,sizeof(double),windowSize,inputFile))==windowSize){   

    for(index1 =0; index1 < WINDOW_SIZE;index1++){
      inputData[index1]*=w[index1];
    // printf("index %d \t %lf\n",index1,inputData[index1]);
    }
     fftw_execute_r2r(plan_r2hc, inputData, windowData);
     for(index1 =0; index1 < windowSize;index1++){
          đầu raData[index1]+=windowData[index1];
     }
     if(chồng chéo!=0)
      fseek(inputFile,(-overlaping)*sizeof(double),SEEK_CUR);
  }
  nếu(i!=0){
     tôi = -i;
    fseek(inputFile,i*sizeof(double),SEEK_END);
    fread(inputData,sizeof(double),-i,inputFile);
    for(index1 =0; index1 < -i;index1++){
      inputData[index1]*=(1.0 - cos(2.0 * M_PI * index1/(windowSize-1)) * 0.5);
    // printf("index %d \t %lf\n",index1,inputData[index1]);
    }
     fftw_execute_r2r(plan_r2hc, inputData, windowData);
     for(index1 =0; index1 < windowSize;index1++){
          đầu raData[index1]+=windowData[index1];
     }

  }


  powVal = outData[0]*outputData[0];
  powVal /= (windowSize*windowSize)/2;
  chỉ số1 = 0;
  fprintf(outputFile,"%lf ",powVal);
  printf(" PSD \t %lf\n",powVal);
  cho (index1 =1; index1<=windowSize;index1++){
            powVal = outData[index1]*outputData[index1]+outputData[windowSize-index1]*outputData[windowSize-index1];
            powVal/=(windowSize*windowSize)/2;
          // powVal = 20*log10(fabs(powVal));
            fprintf(outputFile,"%lf ",powVal);
            printf(" PsD %d \t %10.5lf\n",index1,powVal);
    }


  fftw_free(inputData);
  fftw_free(outputData);
  fftw_free(windowData);
  fclose(inputFile);
  fclose(outputFile);
}

Tại sao kết quả là tôi chỉ nhận được 0,0000? Tôi đã hỏi rất nhiều câu hỏi về windows và cách sử dụng nó, nhưng tôi thực sự không hiểu mình đang làm gì sai ở đây, có ý kiến ​​gì không?

2.Cập nhậtTheo câu trả lời của SleuthEye, kết quả có vẻ tốt, hãy so sánh kết quả với những gì MATLAB đã sử dụng:

  [đầu ra,f] = pwelch(đầu vào,hann(8192));
    cốt truyện (đầu ra);

Nhập khẩu c Sau kết quả của chương trình, PSD giống nhau nhưng tỷ lệ khác nhau:

.

Như bạn có thể thấy, tỷ lệ là khác nhau.

câu trả lời hay nhất

Như chux đã đề cập,dữ liệu đầu ra Các phần tử của mảng không được khởi tạo.

此外，使用 Phương pháp Bartlett Ước tính mật độ phổ công suất thu được sẽ lấy trung bình công suất của các giá trị phổ (thay vì tính công suất trung bình):

outData[index1] += windowData[index1]*windowData[index1];

Cuối cùng, nếu việc chia tỷ lệ phổ là quan trọng đối với ứng dụng của bạn (tức là nếu bạn cần nhiều hơn cường độ tương đối của các thành phần tần số), thì bạn cũng nên đưa hiệu ứng cửa sổ vào hệ số chuẩn hóa, ví dụ: Công thức số được mô tả trong :

Wss gấp đôi = 0,0;
cho (i=0; i
  Wss += w[i]*w[i];
}
Wss *= WINDOW_SIZE;

Vì vậy, hãy tổng hợp tất cả lại và xem xét Đóng gói nửa phức tạpĐịnh dạng bạn đã sử dụng:

đầu raData = fftw_malloc(sizeof(double)*(windowSize/2 + 1));
for (index1 =0; index1 <= windowSize/2; index1++) {
  dữ liệu đầu ra [chỉ mục1] = 0,0;
}
Wss = 0,0;
cho (i=0; i
  Wss += w[i]*w[i];
}
Wss *= WINDOW_SIZE;
...

đếm = 0;
while((i=fread(inputData,sizeof(double),windowSize,inputFile))==windowSize) {   
  for(index1 =0; index1 < WINDOW_SIZE;index1++){
    inputData[index1]*=w[index1];
  }
  fftw_execute_r2r(plan_r2hc, inputData, windowData);

  đầu raData[0] += windowData[0]*windowData[0];
  for(index1 =1; index1 < windowSize/2; index1++)
  {
    double re = windowData[index1];
    double im = windowData[windowSize-index1];
    outData[index1] += re*re + im*im;
  }
  outData[windowSize/2] += windowData[windowSize/2]*windowData[windowSize/2];
  count++;
}
...
nếu (halfSpectrum){
  định mức = số lượng*Wss/2;
  powVal = đầu raData[0]/(2*norm);
  fprintf(outputFile,"%lf ",powVal);
  cho (index1 =1; index1
    powVal = dữ liệu đầu ra[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
  powVal = outData[windowSize/2]/(2*norm);
  fprintf(outputFile,"%lf ",powVal);
}
khác{
  định mức = số lượng*Wss;
  cho (index1 =0; index1<=windowSize/2;index1++){
    powVal = dữ liệu đầu ra[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
  for (index1 =windowSize/2-1; index1>0;index1--){
    powVal = dữ liệu đầu ra[index1]/norm;
    fprintf(outputFile,"%lf ",powVal);
  }
}

Cập nhật (giải thích sự khác biệt với Matlab pwelch):

Theo Octave tiếng vỗ tay Trợ giúp (phải khớp với đầu ra của Matlab):

Mật độ phổ là giá trị trung bình của biểu đồ, được chia tỷ lệ sao cho diện tích dưới phổ giống với bình phương trung bình của dữ liệu.

nói cách khác,

Và hệ số tỷ lệ được cung cấp ở trên áp dụng cho định nghĩa trong đó tỷ lệ sao cho tổng các giá trị phổ rời rạc giống như bình phương trung bình của dữ liệu (nhưcâu hỏi trước đâyphù hợp với mong đợi):

Vì vậy, sự khác biệt trong định nghĩa đưa ra một bổ sung WINDOW_SIZE*tốc độ lấy mẫu các yếu tố (lưu ý rằng nếu không được chỉ định tiếng vỗ tay Đối với tham số thứ tư, bạn sẽ sử dụng mặc định tỷ lệ lấy mẫu tức là 1Hz).

Vì vậy, để làm cho đầu ra nửa phổ của phiên bản C phù hợp tiếng vỗ tay, bạn cần:

định mức = số lượng*Wss/(2*WINDOW_SIZE*sampling_rate);
powVal = đầu raData[0]/(2*norm);
fprintf(outputFile,"%lf ",powVal);
cho (index1 =1; index1
  powVal = dữ liệu đầu ra[index1]/norm;
  fprintf(outputFile,"%lf ",powVal);
}
powVal = outData[windowSize/2]/(2*norm);
fprintf(outputFile,"%lf ",powVal);

Nếu không, hãy chia tỷ lệ tiếng vỗ tay Đầu ra khớp với định nghĩa được sử dụng trong chương trình C:

% Đối với sample_rate tùy ý:
%[output,f] = pwelch(input,hann(8192),[],[],sampling_rate)/(8192*sampling_rate);
% giúp đơn giản hóa các điều sau bằng cách đặt samples_rate = 1
[đầu ra,f] = pwelch(đầu vào,hann(8192))/8192;
cốt truyện (đầu ra);

Về c - cần trợ giúp để chương trình này tính toán PSD, chúng tôi đã tìm thấy một câu hỏi tương tự trên Stack Overflow: https://stackoverflow.com/questions/25142117/