Преобразование Фурье: получение магнитного поля + фаза, а затем использование их для построения исходного сигнала

Привет, парень, я работаю над простыми сигналами, и я хочу рассчитать преобразование Фурье сигнала, получить амплитуду и фазы, а затем восстановить исходный сигнал из этого.

Я основываю свой код на этой теме.

Код:

>> n=0:99;
>> N=length(n);
>> x = sin((2*pi/N).*n).*cos((pi/N).*n);
>> F = fft(x);
>> mag =  sqrt(real(F).^2 + imag(F).^2);
>> phase = atan2(imag(F),real(F));
>> re = mag .* cos(phase);
>> im = mag .* sin(phase);
>> F_i = re + 1i*im;
>> x_i = ifft(F_i);
>> figure;stem(x);figure;stem(x_i);

Я получаю совершенно другие графики.

РЕДАКТИРОВАТЬ: на самом деле я делаю это, чтобы проверить, что произойдет с сигналом, если изменится фаза. Таким образом, для повторного построения сигнала мне понадобится фазовый угол.

Я все еще новичок в Фурье + Matlab, поэтому прошу прощения, если я делаю какую-то случайную глупую ошибку. Я буду признателен, если вы, ребята, сможете указать мне правильное направление. Спасибо.


matlab fft ifft
person kir    schedule 13.11.2013    source источник
comment
Мнимая часть F_i в основном является ошибкой округления, если вы используете F_i=re, цифра воспроизводится правильно. Я недостаточно знаю о преобразованиях Фурье, чтобы сказать вам, почему это так склонно к численной нестабильности.   -  person David    schedule 13.11.2013
comment
Это правда. На самом деле я делаю это для эксперимента, где я изменяю фазовый угол, реконструирую сигнал и смотрю, что изменится. Есть ли другой способ сделать это тогда?   -  person kir    schedule 13.11.2013

Ответы (2)


arrow_upward
2
arrow_downward

Проблема вызвана ошибками округления в phase, поэтому не используйте их при вычислении синуса и конуса фазовых углов. Вместо этого используйте идентификаторы триггеров cos(atan(A))=(1+A^2)^(-1/2) и sin(atan(A))=A*(1+A^2)^(-1/2), и т.д.

re = mag .* real(F)./sqrt(real(F).^2+imag(F).^2);
im = mag .* imag(F)./sqrt(real(F).^2+imag(F).^2);

РЕДАКТИРОВАТЬ: я думаю, что если вы хотите изменить фазовый угол на S, это поможет:

re = mag .* (real(F)*cos(S)-imag(F)*sin(S))./sqrt(real(F).^2+imag(F).^2);
im = mag .* (real(F)*sin(S)+imag(F)*cos(S))./sqrt(real(F).^2+imag(F).^2);

EDIT2: иногда вы все равно будете получать плохие результаты с ненулевой мнимой частью (например, если S=pi), и вам нужно будет построить либо stem(real(x_i)), либо stem(1:length(x_i),x_i), как предложил Луис.

person David    schedule 13.11.2013
comment
Спасибо, это правильно воспроизводит сигнал. Если я хочу изменить im, скажем, добавив константу 10, нужно ли мне делать какие-либо вычисления или будет достаточно im = im + 10? - person kir; 13.11.2013
comment
Да, если вы хотите изменить im, просто сделайте это. Я отредактировал свой ответ, включив в него формулы для re и im с фазовым сдвигом S радиан. - person David; 13.11.2013
comment
@David Предположим, если я хочу изменить всю фазовую матрицу на что-то, скажем, фазово-новое, а затем хочу применить ifft () с использованием mag и фазово-новый, то как мне поступить? - person Jerky; 16.11.2014
comment
@Jerky Возможно, лучше задать новый вопрос, но я не думаю, что мой метод сработает, так как это просто постоянное изменение фазы. Просто замените рассчитанную phase другой матрицей, и она должна работать. - person David; 17.11.2014

arrow_upward
2
arrow_downward

В БПФ такой численной нестабильности нет. Проблема в том, что ошибки округления дают очень маленькую мнимую часть восстановленного сигнала x_i. Это нормально. Действительная часть x_i правильно воспроизводит x, тогда как мнимая часть x_i очень мала. Вы можете проверить с stem(real(x_i)) и stem(imag(x_i))

Теперь stem(x_i) с комплексом x_i отображает мнимую часть в сравнении с реальной частью. С другой стороны, stem(1:length(x_i),x_i) (или, конечно, stem(real(x_i))) отображает действительную часть x_i. Обратите внимание, что это согласуется с поведением plot (см. также ответ на этот вопрос )

Так что просто нарисуйте реальные мнимые части x_i отдельно с stem. Вы увидите, что действительная часть воспроизводит x, как и ожидалось, а мнимая часть пренебрежимо мала (порядка eps).

person Luis Mendo    schedule 13.11.2013
comment
Да, это правда. Интересно, почему мое решение решает проблему? - person David; 13.11.2013
comment
@David Каким-то образом это не вызывает ошибок округления, следовательно, нет мнимой части в stem и нет проблем с построением графика. Странная часть заключается в том, почему stem делает это со сложным аргументом. Было бы больше смысла, если бы он игнорировал мнимую часть, как это делает plot(x,y) - person Luis Mendo; 13.11.2013
comment
@David На самом деле stem(x,y) действительно ведет себя как plot(x,y): игнорирует мнимую часть y. - person Luis Mendo; 13.11.2013