Matlab: коэффициенты фурье функции sign () колеблются

Я пишу функцию, которая принимает изменяющийся во времени сигнал и возвращает БПФ. Это следует за документацией Matlab:

Он хорошо работает для действительно простых сумм синусоид.

Для других смоделированных сигналов, таких как пошаговые функции, я получаю что-то fubar, которое колеблется в виде пилообразной формы. Вот минимальный пример с использованием ступенчатой ​​функции тяжелой стороны, которая отображает мнимую составляющую веса БПФ:

%Create and plot sign function 
Fs = 1000;
dt = 1/Fs;
tvals = -1: dt: 1-dt;
yvals = sign(tvals); %don't use heaviside() it is very slow

subplot(2,1,1);
    plot(tvals, yvals, 'k', 'LineWidth', 2); hold on; 
    axis([-1 1 -1.1 1.1]);

%Calculate center-shifted FFT and plot imaginary part
fftInitial = fft(yvals);
n = length(fftInitial);
dF = Fs/n;
frequenciesShifted = -Fs/2: dF: Fs/2-dF;  %zero-centered frequency range
fftShifted = fftshift(fftInitial/length(yvals));

subplot(2,1,2);
plot(frequenciesShifted, imag(fftShifted), 'b', 'LineWidth', 2);hold on
xlim([-8 8])

Вот итоговый сюжет: enter image description here

Обратите внимание, что воображаемое решение известно как 2/jw или j(-2/w).

Обратите внимание, помимо пилообразного зуба, который является моей основной проблемой, конверт этих весов, похоже, даже не следует за этим. На самом деле это кажется перевернутым с точки зрения происхождения. Не уверен, какое сочетание вещей я сделал здесь неправильно.

Основываясь на некоторых полезных отзывах, люди указали на этот вопрос:
Аналитическое преобразование Фурье против БПФ функций в Matlab

В частности, отсутствие нуля во временном массиве могло вызвать проблемы, что было главной проблемой. Массив времени в моем коде включает ноль, поэтому он не похож на дубликат. Я переместил свой ноль в начало массива, сдвигая шаг по времени (хотя, честно говоря, я сделал много fft (), не делая этого, и все они выглядели нормально, поэтому я не думаю, что это проблема, но Я просто хочу пока убрать это как проблему). Итак, мы получаем:

Fs = 1000;
dt = 1/Fs;
tvals = 0: dt: 2-dt;
yvals = sign(tvals-1); %don't use heaviside() it is very slow
zeroInd = find(yvals == 1, 1, 'first');
yvals(zeroInd-1) =0;
%Calculate center-shifted FFT and plot imaginary part
fftInitial = fft(yvals);
n = length(fftInitial);
dF = Fs/n;
frequenciesShifted = -Fs/2: dF: Fs/2-dF;  %zero-centered frequency range
fftShifted = fftshift(fftInitial/length(yvals));

%Plot stuff
subplot(2,1,1);
plot(tvals, yvals, 'k', 'LineWidth', 2); hold on; 
axis([0 2 -1.1 1.1]);  
subplot(2,1,2);
plot(frequenciesShifted, imag(fftShifted), 'b', 'LineWidth', 2);hold on
xlim([-8 8])
grid on;

Я по-прежнему получаю ту же зигзагообразную функцию с тем же неправильным конвертом. Итак, хотя я вижу, что мой вопрос и этот вопрос тесно связаны, я не уверен, что они дублируются. И я действительно хотел бы иметь возможность нарисовать наполовину разумный график значений этих компонентов (я работаю с некоторыми физиологическими сигналами, которые в основном являются ступенчатыми функциями (они перемещают очень быстрее по сравнению с моими измерительными приборами, так что это не просто академическое упражнение для меня).

Повторяющийся вопрос использует другую функцию, но имеет точно такую ​​же проблему, с которой вы столкнулись: БПФ предполагает, что начало координат (t = 0) является первым элементом во входном массиве. Вы вычисляете БПФ сдвинутой ступенчатой ​​функции.

Cris Luengo 13.09.2018 19:46

Этот вопрос на самом деле не является дубликатом: во-первых, у меня есть ноль в моем представлении времени, и это было проблемой в их сообщении. Во-вторых, моя основная проблема - это колебательное поведение FT, а не зеркальное отражение, как указано в моем заголовке. Но я отредактирую свой вопрос, чтобы упомянуть этот пост.

eric 13.09.2018 20:06

@CrisLuengo Cris Я прочитал это: ... Как вы там написали, «БПФ ожидает, что исходная точка будет в первом (крайнем левом) образце. Это то, для чего предназначен ifftshift:« Итак, то, что я сделал в моем исправленном вопросе, было просто сделать 0 мой первый момент времени. Я думал, что это избавит от необходимости использовать iffshift. Теперь посмотрим на это более внимательно.

eric 13.09.2018 20:48

@CrisLuengo, почему это чрезвычайно простое решение работает для функции единичного импульса (без необходимости возиться со временем и т. д.) С великолепным сюжетом, но когда я делаю то же самое с одношаговой функцией, это катастрофа? stackoverflow.com/questions/19717779/… То есть почему там ответ такой простой и легкий. Я просто не вижу разницы (особенно учитывая, что импульс - это просто сумма двух сигналов!).

eric 13.09.2018 21:25

Я был неправ, у вас есть эта проблема, но есть и вторая. Извините. Смотрите мой ответ. Вопрос с единичным импульсом, который вы связываете здесь, отображает величину преобразования Фурье, поэтому сдвиги в сигнале не видны.

Cris Luengo 13.09.2018 21:29
1
5
510
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

В вашем коде есть две проблемы:

  1. ДПФ (алгоритм БПФ вычисляет ДПФ) берет начало координат в крайнем левом интервале. Создание yvals таким образом, что начало координат находится посередине, приводит к тому, что выходной частотный спектр становится спектром сигнала, сдвинутого на половину его длины. Это приводит к колебаниям очень высокой частоты. Исправление состоит в том, чтобы использовать ifftshift для входных данных до, вызывая fft. Посмотреть еще в этом другом вопросе

  2. ДПФ предполагает (можно интерпретировать как допущение), что входной сигнал является периодическим. Это приводит ко второму большому прыжку. По сути, ваш сигнал выглядит как функция сдвинутого блока, и, следовательно, ваше преобразование выглядит как функция sinc с измененной фазой. Решение состоит в том, чтобы применить оконную функцию к вашему входу до, вызывая fft. См., Например, этот другой вопрос.

Измените свой код следующим образом:

yvals = sign(tvals);
yvals = yvals .* hanning(numel(yvals), 'periodic').'; % Apply windowing function

% ...

fftInitial = fft(ifftshift(yvals)); % Shift signal before calling FFT

Вот результат, который теперь дает ваш код:

graph

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

eric 13.09.2018 21:59

Другие вопросы по теме