Использование БПФ и правильное масштабирование по Гауссу в MATLAB - PullRequest
Новая
       75

Использование БПФ и правильное масштабирование по Гауссу в MATLAB

0 голосов
/ 10 июля 2020

У меня проблема с получением правильного вывода для преобразования Фурье гауссова exp (-2 (x ^ 2 + y ^ 2). Я знаю аналитически, что преобразование Фурье должно быть 0,25 * exp (-0,125 (k ^ 2) + j ^ 2)), где k и j - переменные Фурье, но выход абсолютного значения на выходе БПФ не соответствует этому. Это решение аналити c:

по сравнению с моим выводом:

Это соответствующий код: 

Fs = 10;                                        %space frequency
Dx = 1 / Fs;                                    %sampling period
x = -10: Dx : 10;                               %space vector on domain [-10, 10], one dimensional
L = length(x);                                  %length of signal
[X, Y] = meshgrid(x, x);                        %2D domain defined by the space vector x

U =  exp(-2 *( X.^2 + Y.^2));                                      
n = 2^nextpow2(L);                              %padding for FFT to optimize performance

U_hat = fft2(U, n, n);                          %FFT of the Gaussian

Dk = Fs*((-n/2):((n/2) - 1)) / n;               %space frequency domain (ie, the fourier domain) in one dimension, shifted, rescaled by n
P = abs(fftshift(U_hat) / n);                   %power spectrum (ie | X_hat | = sqrt(X_hat * complex_conjugate(X_hat)), shifted
                                                                                                                 
[K1, K2] = meshgrid(Dk, Dk);                    %Fourier domain

mesh(K1, K2, P);
title('Gaussian Pulse in Frequency Domain');
xlabel('Frequency (k_1)');
ylabel('Frequency (k_2)');
zlabel('|P(f)|');

Я испортил частотную область?

Ответы [ 2 ]

1 голос
/ 11 июля 2020

Благодарим Криса Луена go за ответ в комментариях - решение analyti c, представленное в вопросе, было неверным. Mathematica сообщила, что преобразование Фурье - 0.25 * exp(-0.125(k^2 + j^2)), но фактическое преобразование - (pi / 2) * exp(-( (pi^2) / 2) * (k^2 + j^2)). Код правильный как есть.

0 голосов
/ 11 июля 2020

Я довольно долго анализировал вашу реализацию. Ваша реализация в частотной области выглядит хорошо. Почему?

  1. У вас есть гауссовский как FFT O / P для гауссовского как ввода. (Подробнее об этом позже)
  2. Приращение элемента частотного спектра, т.е. df = Fs / n. Это говорит о том, что по оси X вашего частотного преобразования нет ничего неправильного.
  3. Ваша ось Y преобразования равна абс. Fftshift или fft, деленному на n, чтобы нормализовать преобразование. Вы не умножили на 2, что большинство людей приняли бы, поскольку вас интересуют также отрицательные частоты.

Итак, реализация верна, но конечный результат неверен ... Это может означать только один вещь. Ваш гауссов сигнал выглядит только гауссовским, но на самом деле это не так! Если ваши средние значения равны 0, а дисперсия и sd равны 1. Оно просто становится

U = (1 / (2 pi)) exp (- (X. ^ 2 + Y. ^ 2) / 2) --- (1)

Даже если вы проигнорируете первую часть (1), то есть 1/2 * pi как есть простое масштабирование амплитуды, ваше уравнение * exp (-2 (X. ^ 2 + Y. ^ 2)) все еще не соответствует основному c уравнению гауссовой волны .

Пожалуйста, еще раз проверьте введенные данные.

...