2015-05-22 3 views
2

Я пишу простой код для сложного сложного DFT в c с библиотекой fftw3. Я написал файл с двойными данными входного массива, поэтому я могу сравнить с функцией matlab fft. Я пытаюсь выполнить обратное преобразование из массива преобразования, но результаты и первый входной массив различны. это мои результаты:fftw3 обратное преобразование не работает

FFTW3 TRANSFORM WELCOME <<<<< 

enter the number (integer) N of samples (Bit: 64) (preferably power of 2):8 

SAMPLE INPUT 
in[0][0] = -216448918.015237  in[0][1] = 0.000000 
in[1][0] = 948904790.062151   in[1][1] = 0.000000 
in[2][0] = 826811206.185300   in[2][1] = 0.000000 
in[3][0] = 1868763250.342451  in[3][1] = 0.000000 
in[4][0] = 703135606.077152   in[4][1] = 0.000000 
in[5][0] = -1989016445.622210  in[5][1] = 0.000000 
in[6][0] = 1912963650.704585  in[6][1] = 0.000000 
in[7][0] = 811527262.805480   in[7][1] = 0.000000 

Hit enter to continue ... 


FORWARD TRANSFORM COEFFICIENTS 

out[0][0] = 4866640402.539672  out[0][1] = 0.000000 
out[1][0] = 410260768.150135  out[1][1] = -1738850319.926936 
out[2][0] = -2253088168.827970  out[2][1] = 3720402168.707990 
out[3][0] = -2249429816.334913  out[3][1] = -3911155208.965507 
out[4][0] = 1586282687.363928  out[4][1] = 0.000000 
out[5][0] = -2249429816.334913  out[5][1] = 3911155208.965507 
out[6][0] = -2253088168.827970  out[6][1] = -3720402168.707990 
out[7][0] = 410260768.150135  out[7][1] = 1738850319.926936 
do you want to calculate the inverse-transform? (y/n) 
y 


INVERSE TRANSFORM COEFFICIENTS 
rev[0][0] = -1731591344.121896  rev[0][1] = 0.000000 
rev[1][0] = 7591238320.497208  rev[1][1] = 0.000000 
rev[2][0] = 6614489649.482399  rev[2][1] = 0.000000 
rev[3][0] = 14950106002.739609  rev[3][1] = 0.000000 
rev[4][0] = 5625084848.617215  rev[4][1] = 0.000000 
rev[5][0] = -15912131564.977680  rev[5][1] = 0.000000 
rev[6][0] = 15303709205.636681  rev[6][1] = 0.000000 
rev[7][0] = 6492218102.443840  rev[7][1] = 0.000000 

Как вы видите, «в» и «REV» массивы различны, но прямое преобразование является правильным. Я сравнил его с Matlab, и результаты остались прежними. Когда я выполняю обратное преобразование с помощью matlab, я получаю входной массив. Что мне делать?

это мой с код:

#include <fftw3.h> 
#include <math.h> 
#include <stdio.h> 
#include <stdlib.h> 
#define PI 3.141592653589 

int main() 
{ 

    fftw_complex *in, *out, *rev; 
    int i,f0,A,N; 
    char no; 
    FILE *fp; 
    fftw_plan p; 

    printf("\n\n>>>>> FFTW3 TRANSFORM WELCOME <<<<<"); 
    printf("\n\n enter the number (integer) N of samples (bit: %ld) (preferably power of 2):",(sizeof(fftw_complex)/2)*8); 
    scanf("%d",&N); 

    //f0 = 50; 
    //A = 1; 


    //allocating memory for input & output arrays 
    in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N); 
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N); 
    rev = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N); 



    //Opening the data file 
    if((fp=fopen("lista_numeri_double.dat","rb"))==NULL) 
    { 
    printf("\nError reading file\n"); 
    exit(1); 
    } 


    printf("\nSAMPLE INPUT"); 
    //assigning samples read from the file 
    for(i=0;i<N;i++) 
    { 
    //in[i][0] = A * cos(2*PI*f0*i/N); 
    fread(&in[i][0],sizeof(double),1,fp); 
    in[i][1]=0; 

    printf("\nin[%d][0] = %f \t\tin[%d][1] = %f",i,in[i][0],i,in[i][1]); 
    } 


    //plan and execute transform 
    p = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_ESTIMATE); 
    fftw_execute(p); 

    printf("\n\n Hit enter to continue ... \n"); 
    scanf("%c",&no); 


    //print output values 
    printf("\n\nFORWARD TRANSFORM COEFFICIENTS\n"); 
    for(i=0;i<N;i++) 
    { 
    printf("\nout[%d][0] = %f \t\tout[%d][1] = %f",i,out[i][0],i,out[i][1]); 
    } 

    fftw_destroy_plan(p); 

    printf("\n do you want to calculate the inverse-transform? (y/n) \n"); 
    scanf ("%c",&no); 

    if(no=='y') 
    { 

    //plan and execute inverse transform 
    p = fftw_plan_dft_1d(N,out,rev,FFTW_BACKWARD,FFTW_ESTIMATE); 
    fftw_execute(p); 

    printf("\n\nINVERSE TRANSFORM COEFFICIENTS\n"); 
    for(i=0;i<N;i++) 
    { 
     printf("rev[%d][0] = %f \t\trev[%d][1] = %f\n",i,rev[i][0],i,rev[i][1]); 
    } 

    fftw_destroy_plan(p);   
    } 

    return 0; 
} 

ответ

2

Разница между Matlab и FFTW поставляется с коэффициентом масштабирования применяется для преобразования.

Принимая во внимание, что Matlab's FFT нормирован, алгоритм, используемый FFTW, как описано in FFTW's documentation, не нормирован. Другими словами, преобразование полного круга с использованием FFTW (вперед с последующим обратным) масштабирует результат в N.

Соответственно, сравнивая массив in и rev показывает, что rev масштабируются последовательными 8 раз (размером N из преобразования, используемых в вашем примере).

+0

Я думал, что в эти дни matlab украл 'fftw' и включил его в свои собственные. [3] FFTW (http://www.fftw.org) [4] Frigo, M. and SG Johnson, «FFTW: адаптивная архитектура программного обеспечения для БПФ», Материалы Международной конференции по акустике, речи, и обработка сигналов, т. 3, 1998, стр. 1381-1384. –

+0

@GAlexander afaik matlab использует модифицированную версию fftw, которая удобно масштабируется для обратной совместимости. – SleuthEye

+0

в pld days fftw будет раздавить matlab fft настолько, что никто, кого я знал, не использовал реализацию matlabs. –

Смежные вопросы