БПФ вещественного симметричного вектора не является вещественным и симметричным.

Question

Мне трудно понять, какой должна быть простая концепция.Я построил в MATLAB действительный и симметричный вектор.Когда я беру БПФ в MATLAB, результат имеет значительную мнимую составляющую, хотя правила симметрии преобразования Фурье гласят, что ПФ действительной симметричной функции также должно быть вещественным и симметричным.Мой пример кода:

N = 1 + 2^8;
k = linspace(-1,1,N);

V = exp(-abs(k));

Vf1 = fft(fftshift(V));
Vf2 = fft(ifftshift(V));
Vf3 = ifft(fftshift(V));
Vf4 = ifft(ifftshift(V));
Vf5 = fft(V);
Vf6 = ifft(V);

disp([isreal(Vf1) isreal(Vf2) isreal(Vf3) isreal(Vf4) isreal(Vf5) isreal(Vf6)])

Результат:

0 0 0 0 0 0

Никаких комбинаций (i)fft или (i)fftshift в результате получим действительный симметричный вектор.Я пробовал как с четным, так и с нечетным N (N = 2^8 против. N = 1+2^8).

я попробовал посмотреть k+flip(k) и имеются остатки порядка eps(1), но остатки также симметричны, и мнимая часть БПФ не выглядит нечеткой порядка eps(1), а скорее с величиной, сравнимой с реальной частью.

Какую ослепительно очевидную вещь я упускаю?

Ослепительно очевидная вещь, которую мне не хватало:

БПФ не является интегралом по все пространство, поэтому предполагается периодический сигнал.Выше я дублирую последний пункт в периоде, когда выбираю четный N, и поэтому нет возможности сдвинуть его, чтобы поставить нулевую частоту в начало без дробной индексации, которой не существует.

Несколько слов о моем выборе k.Это не произвольно.Реальная проблема, которую я пытаюсь решить, — это создать модель интерферограммы FTIR, которую я затем проведу БПФ, чтобы получить спектр. k это расстояние, которое проходит интерферометр, которое преобразуется в частоту в волновых числах.В реальной задаче будут использоваться различные коэффициенты масштабирования, чтобы производящая функция V давала физически значимые числа.

Answer 1

@Ивон абсолютно прав в своем комментарии о симметрии.Ваш входной сигнал выглядит симметричным, но это не потому, что симметрия связана с началом координат 0.С использованием Линспейс в Matlab для построения сигналов — по большей части плохой выбор.Попытка исправить результаты с помощью fftshift также является плохой идеей.

Вместо этого используйте:

k = 2*(0:N-1)/N - 1;

и вы получите ожидаемый результат.Однако мнимая часть преобразованных значений не будет равна нулю.Есть некоторый числовой шум.

>> max(abs(imag(Vf5)))
ans =
2.5535e-15

Ответ на вопрос Ивона:

Почему?>> N = 1+2^4 N = 17 >> x=linspace(-1,1,N) x = -1,0000 -0,8750 -0,7500 -0,6250 -0,5000 -0,3750 -0,2500 -0,1250 0 0,1250 0,2500 0,3750 0,5000 0,6250 0 .7500 0,8750 1,0000 >> y=2*(0:N-1)/N-1 y = -1,0000 -0,8824 -0,7647 -0,6471 -0,5294 -0,4118 -0,2941 -0,1765 -0,0588 0,0588 0,1765 0,2941 0,411 8 0,5294 0,6471 0,7647 0,8824 – Ивон 1

Ваш пример - это не симметричная (четная) функция, а антисимметричная (нечетная) функция.Однако это не имеет никакого значения.

Для антисимметричной функции длины N справедливо следующее утверждение:

f[i] == -f[-i] == -f[N-i]

Индекс i работает от 0 до N-1.

Давайте посмотрим, что произошло с i=2.Помните, счет начинается с 0 и заканчивается 16.

x[2] = -0.75
-x[N-2] == -x[17-2] == -x[15] = (-1) 0.875 = -0.875
x[2] != -x[N-2]

y[2] = -0.7647
-y[N-2] == -y[15] = (-1) 0.7647
y[2] == y[N-2]

Проблема в том, что начало векторов Matlab начинается с 1.Векторы по модулю (периодические) начинаются с начала координат 0.Эта разница приводит ко многим недоразумениям.

Другой способ объяснить, почему linspace(-1,+1,N) неверен:

Представьте, что у вас есть вектор, который имеет один период периодической функции, например, функция Cosinus.Этот единственный период является одним из бесконечного числа периодов.Первое значение вашего вектора косинуса не должно совпадать с последним значением вашего вектора.Однако именно это и делает linspace(-1,+1,N).В результате получается последовательность, в которой последнее значение периода 1 совпадает со значением первой выборки следующего периода 2.Это не то, чего вы хотите.Чтобы избежать этой ошибки, используйте t = 2*(0:N-1)/N - 1.Расстояние t[i+1]-t[i] равно 2/N, и последнее значение должно быть t[N-1] = 1 – 2/N, а не 1.

Ответ на второй комментарий Ивона

Что бы вы ни вводили во входной вектор ДПФ/БПФ, теоретически оно интерпретируется как периодическая функция.Но дело не в этом.

ДПФ выполняет интегрирование.

fft(m) = Sum_(k=0)^(N-1) (x(k) exp(-i 2 pi m k/N )

Первое значение x(k=0) описывает амплитуду первого интервала интегрирования длиной 1/N.Второе значение x(k=1) описывает амплитуду второго интервала интегрирования длиной 1/N.И так далее.

Самый последний интервал интегрирования симметричной функции заканчивается с тем же значением, что и первый образец.Это означает, что отправной точкой последнего интервала интегрирования является k=N-1 = 1-1/N.Ваш входной вектор содержит отправные точки интервалов интегрирования.

Следовательно, последняя точка симметрии k=N является точкой функции, но не является начальной точкой интервала интегрирования и поэтому не является членом входного вектора.

Answer 2

Его

Vf = fftshift(fft(ifftshift(V)));

То есть вам нужно ifftshift во временной области, так что выборки интерпретируются как образцы симметричной функции, а затем fftshift в частотной области, чтобы снова сделать симметрию очевидной.

Это работает только для N странный.Для N даже понятие симметричной функции не имеет смысла:нет способа сдвинуть сигнал так, чтобы он был симметричен относительно начала координат (начало координат должно находиться «между двумя выборками», что невозможно).

Для вашего примера V, приведенный выше код дает Vf действительный и симметричный.Следующий рисунок был создан с помощью semilogy(Vf), так что можно увидеть как маленькие, так и большие значения.(Конечно, вы можете изменить горизонтальную ось так, чтобы график располагался по центру частоты 0, как и должно быть;но в любом случае граф выглядит симметричным.)

Answer 3

У вас возникла проблема при реализации понятия «симметрия».Чисто вещественная, четная (или «симметричная») функция имеет функцию преобразования Фурье, которая также является вещественной и четной.«Четный» — это симметрия относительно оси Y или линии t=0.

Однако при реализации сигнала в Matlab вы всегда начинаете с t=0.То есть невозможно «определить» сигнал, начиная с начала времени.

Недолгий поиск в Интернете привел меня к следующему:Правильное использование fftshift и ifftshift при вводе в fft и ifft..

Как это сделал Луис указал, вам нужно выполнить ifftshift перед подачей сигнала в fft.Причина никогда не была документирована в Matlab, а только в этой теме.По историческим причинам выходные данные AND входы из fft и ifft поменяны местами.То есть вместо заказанного у -N/2 к N/2-1 (естественный порядок), сигнал во временной или частотной области упорядочивается из 0 к N/2-1 а потом -N/2 к -1.Это означает, что правильный способ кодирования fft( ifftshift(V) ), но большинство людей в большинстве случаев игнорируют это.Почему его молча игнорируют, а не создают огромные проблемы, так это то, что больше всего беспокойства вызывает амплитуда сигнала, а не фаза.Поскольку круговое смещение не влияет на амплитудный спектр, это не проблема (даже для ребят из Matlab, написавших документацию).

Чтобы проверить равенство амплитуд -

Vf2 = fft(ifftshift(V));
Vf5 = fft(V);
Va2 = abs(fftshift(Vf2));
Va5 = abs(fftshift(Vf5));

>> min(abs(Va2-Va5)<1e-10)

ans =

     1

Чтобы увидеть, как сильно не в фазе -

Vp2 = angle(fftshift(Vf2));
Vp5 = angle(fftshift(Vf5));

В любом случае, как я писал в комментарии, после копирования и вставки вашего кода в свежий и чистый Matlab, он дает 0 1 0 1 0 0.

Что касается вашего вопроса о N = четном и N = нечетном, я считаю, что когда N = четный, сигнал не симметричен, поскольку по обе стороны от начала координат времени имеется неодинаковое количество точек.

Answer 4

Просто добавьте следующую строку после «k = linspace(-1,1,N);»

k(end)=[];

он удалит последний элемент массива.Это определено как симметричный массив.

также учтите, что isreal(complex(1,0)) является ложным!!!Функция isreal просто проверяет формат хранения памяти.поэтому 1+0i в приведенном выше примере не является реальным.

Вы определили свою функцию для проверки действительных чисел (вот так)

myisreal=@(x) all((abs(imag(x))<1e-6*abs(real(x)))|(abs(x)<1e-8));

Наконец, ваш исходный код должен выглядеть примерно так:

N = 1 + 2^8;
k = linspace(-1,1,N);

k(end)=[];

V = exp(-abs(k));


Vf1 = fft(fftshift(V));

Vf2 = fft(ifftshift(V));
Vf3 = ifft(fftshift(V));
Vf4 = ifft(ifftshift(V));
Vf5 = fft(V);
Vf6 = ifft(V);

myisreal=@(x) all((abs(imag(x))<1e-6*abs(real(x)))|(abs(x)<1e-8));

disp([myisreal(Vf1) myisreal(Vf2) myisreal(Vf3) myisreal(Vf4) myisreal(Vf5) myisreal(Vf6)]);