Результаты FFTW отличаются от FFT в MATLAB

Я сравниваю прямое БПФ, используя FFTW и MATLAB fft. Входной сигнал является гауссовым. Код:

FFTW с использованием C:

float *signal; /*input signal*/
int nt; /*length of the signal*/
fftw_complex *in, *out;
fftw_plan plan1;

in = fftw_malloc(nt*sizeof(fftw_complex)); 
out = fftw_malloc(nt*sizeof(fftw_complex));
for (j=0;j<nt;j++){
        in[j][0]=(double)signal[j];
        in[j][1]=0.0;
}    
plan1 = fftw_plan_dft_1d(nt, in, out, -1, FFTW_ESTIMATE);
fftw_execute(plan1);        
fftw_destroy_plan(plan1);

for (j=0;j<nt;j++){
        real[j]=(float)out[j][0];
        imag[j]=(float)out[j][1];
}

fft функция в MATLAB:

fft(signal);

Я рисую действительную и мнимую части обоих результатов:

fft против fftw

Действительная часть имеет почти одно и то же значение, а мнимая часть имеет совсем другие значения. Как решить эту проблему?


c matlab fft fftw
person clduan    schedule 24.06.2019    source источник
comment
Различия в Imag находятся в диапазоне от 10 до -15 степени. Так что не так уж и отличается.   -  person Kami Kaze    schedule 24.06.2019
comment
Другой факт, потому что мой сигнал - это действительные числа, но в FFTW я создаю сложное ДПФ, что означает, что сигнал сложный с нулевыми мнимыми значениями. Для реального сигнала out[i] является сопряженным out[n-i], как показывает результат Matlab fft. Но образ FFTW так себя не ведет.   -  person clduan    schedule 24.06.2019
comment
Хорошим тестом было бы применение обратного преобразования в обоих случаях и сравнение результатов, это покажет вам, имеют ли значение различия в мнимой части (я также думаю, что это не так).   -  person schnaader    schedule 24.06.2019
comment
что означает, что сигнал является комплексным с нулевыми мнимыми значениями. ДПФ имеет нулевые мнимые значения, потому что входные данные симметричны, а не потому, что входные данные действительны. Кроме того, зная, что мнимые значения должны быть равны нулю, вы должны были заметить, что они очень близки к нулю. Вы просто смотрите на ошибки округления с плавающей запятой здесь. Ошибки округления не обязательно должны быть симметричными.   -  person Cris Luengo    schedule 24.06.2019
comment
Я также вижу, что вы неправильно создали свой сигнал Гаусса, поэтому ваш ДПФ не похож на сигнал Гаусса. Пространственная область также имеет начало в первом элементе, вам нужно циклически сдвинуть ваш сигнал так, чтобы центр гауссианы находился в первой выборке, а левая сторона гауссианы приходилась на правый конец сигнала.   -  person Cris Luengo    schedule 24.06.2019
comment
Да, если я сделаю обратное преобразование в обоих случаях, результаты будут одинаковыми. Но если я умножу фильтр на ряд Фурье и обратное преобразование, то будет немного другое из-за не точного сопряжения мнимой части.   -  person clduan    schedule 24.06.2019
comment
Вы можете решить не точную сопряженную часть, обнулив все мнимые значения, меньшие эпсилон (например, 10 ^ -14), поскольку вы все равно знаете, что это шум с плавающей запятой.   -  person schnaader    schedule 24.06.2019
comment
@schnaader Это хороший способ. Спасибо.   -  person clduan    schedule 24.06.2019

Ответы (2)


arrow_upward
2
arrow_downward

Вы должны смотреть на масштабный коэффициент графика слева над графиком «Имаг». Там написано 10^-15. Это довольно мало по сравнению с реальной величиной сигнала (по крайней мере, большая часть, которая составляет> 10 ^ 1), поэтому результаты очень похожи.

Алгоритмы с плавающей запятой в целом имеют тенденцию не давать точно такой же результат, если они не реализованы точно таким же образом. (Да и то они могут отличаться разными вариантами округления).

Этот QA может дать некоторое представление: примеры неточности с плавающей запятой

person Kami Kaze    schedule 24.06.2019
comment
О совсем малом в таком случае полезно судить только относительно чего-то другого. Абсолютной величины недостаточно. В этом случае, однако, любая из нескольких характеристик реального сигнала (максимальный, средний, абсолютный диапазон) дает подходящую основу для сравнения, чтобы сделать вывод, что да, мнимый сигнал в данном случае пренебрежимо мал. - person John Bollinger; 24.06.2019
comment
@JohnBollinger Совершенно верно, я только что сделал вывод о реальной части спектра. Сделал шанс сделать это яснее. - person Kami Kaze; 24.06.2019
comment
Алгоритмы с плавающей запятой, как правило, не дают точно такой же результат, если они не реализованы точно таким же образом. Более того, FFTW может давать разные результаты в зависимости от конкретного решения (FFTW). план) выбранный метод. Например, вычисление ДПФ массива данных из 40 элементов будет несколько отличаться, если он разбит на подэлементы 4x10 или 5x8. Запустите его на разных компьютерах или просто получите другую схему памяти, и вы можете получить разные результаты. Чтобы получить последовательные* результаты от FFTW, необходимо сохранить планы и использовать их повторно. - person Andrew Henle; 24.06.2019
comment
* - последовательно, так как одни и те же входные данные всегда производят одни и те же выходные данные до последнего бита. - person Andrew Henle; 24.06.2019
comment
@AndrewHenle Приятно знать, что я никогда не работал с FFTW, поэтому я не могу касаться этого вопроса, но обычно это подпадает под то, что я сказал о реализации различных алгоритмов, поскольку у FFTW, похоже, есть несколько вариантов на выбор. Этого кто-то может не ожидать такое поведение, так что спасибо за дополнение. - person Kami Kaze; 24.06.2019

arrow_upward
0
arrow_downward

Округлив до ближайших 0,001% полной шкалы (действительной), обратите внимание, что все мнимые значения равны нулю.

person hotpaw2    schedule 24.06.2019