Matlab: conv() -> fft() * fft() -> ifft()

Добрый день всем!

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

Но результаты, которые я получаю, как-то совершенно неверны, и, вероятно, здесь я сочетаю разные неправильные шаги. Я уже просматривал подобные темы здесь и на других сайтах, таких как developmentpez, но не смог понять причину. Буду признателен за любую помощь.

Допустим, мой истинный сигнал f[.] — это просто импульс в момент времени 1, а импульсная характеристика g[.] является гауссовой. Я вычисляю их свертку h[.] по conv(), а затем, по сути, хочу найти ifft( fft[h]./fft[g] ), ожидая, что это будет f[.].

Первая проблема заключается в том, что conv() создает массив из элементов n+m-1, где n,m — длины массивов аргументов. Итак, чтобы выполнить fft[h]./fft[g], мне нужно сделать что-то с длиной g. Это первое подозрительное место, где я могу действовать неправильно (см. код). Как правильно это сделать?

Вторая проблема заключается в том, что я получаю что-то сильно отличающееся от исходного истинного сигнала.

Третья проблема заключается в том, что я не могу понять, как принимать сдвиги сигнала. В матлабе мне приходится оперировать с положительными по времени сигналами, но, например, гауссовская импульсная характеристика имеет как отрицательные по времени, так и положительные по времени элементы, поэтому, чтобы работать с ней здесь, мне нужно сдвинуть ее «вперед» ( взгляд будет двигаться вправо), и чем мне нужно «сдвинуть» результат?

Спасибо!

Вот моя хрень по этому поводу :)

close all;

TrueSignal = zeros( 101, 1 ); % impulse in t = 1.
TrueSignal( 1 ) = 1;
ImpulseResp = normpdf(-1:0.02:1)/normpdf( 0 ); % 101 elements array

figure;
subplot( 2,2,1 );
title('True signal')
plot( TrueSignal );
subplot( 2,2,2 );
title('Impulse response')
plot( ImpulseResp );

Conv = conv( TrueSignal, ImpulseResp ); % produces 201 elements array.
subplot( 2,2,3 );
title('Convolution')
plot( Conv );

% Wrong? I need a 201 elements array to represent the impulse response.
ImpulseResp_sparse = normpdf( -1:0.01:1 )/normpdf( 0 );
FIR = fft( ImpulseResp_sparse )/201;

Inverse = ifft( fft( Conv )./FIR ); % UPD Added fft() according to one of comments, bad mistake, but still not preventing.

subplot( 2,2,4 );
title('What is that???')
plot( abs( Inverse ) ); % It's weird! With no abs(), result is even more weird!

matlab convolution signal-processing fft
person agronskiy    schedule 11.03.2013    source источник
comment
БПФ предполагает, что ваш сигнал является периодическим. Вы можете справиться с этим, используя соответствующие отступы. Это обсуждается в главе о БПФ в Численные рецепты, и здесь она может оказаться полезной.   -  person sfstewman    schedule 12.03.2013
comment
См. также объяснение круговой свертки в Википедии, особенно диаграмму. Использование БПФ/ОБПФ для вычисления сверток возможно, но вы должны обратить внимание на детали. Возможно, вам больше повезет с этим вопросом на dsp.stackexchange.com.   -  person mtrw    schedule 12.03.2013
comment
Связанный вопрос: Проверьте теорему свертки   -  person Eitan T    schedule 12.03.2013
comment
Всем спасибо! Я получил много информации для размышления, которая, безусловно, была очень полезной!   -  person agronskiy    schedule 12.03.2013

Ответы (2)


arrow_upward
2
arrow_downward

Прямое использование fft для свертки приведет к круговой свертке, тогда как то, что вы хотите (и что conv делает) — это линейная свертка. Таким образом, чтобы реализовать такую ​​схему с fft, вам придется дополнить сигналы нулями до длины m+n-1.

Вот пример, показывающий эквивалентность выходных данных линейной свертки на основе conv и fft:

x=rand(4,1);y=rand(3,1); %sample data
out1=conv(x,y);          %output from conv()
X=fft(x,6);Y=fft(y,6);   %zero pad and compute fft
out2=ifft(X.*Y);         %output from fft based lin. conv.

Вы можете проверить, что out1 и out2 совпадают (с точностью до FP).

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

person abcd    schedule 11.03.2013

arrow_upward
2
arrow_downward

  1. Вы можете использовать filter(g, 1, f) вместо conv(g, f), чтобы избежать проблем с результирующей длиной. Или вырезать полученный массив.
  2. Похоже, вы забыли сделать один ффф: Inverse = ifft( fft(Conv)./FIR );

Это работает для меня:

   close all;

   TrueSignal = zeros( 101, 1 ); % impulse in t = 1.
   TrueSignal( 1 ) = 1;
   ImpulseResp = normpdf(-1:0.02:1)/normpdf( 0 ); % 101 elements array

   figure;
   subplot( 2,2,1 );
   title('True signal')
   plot( TrueSignal );
   subplot( 2,2,2 );
   title('Impulse response')
   plot( ImpulseResp );

   Conv = filter( TrueSignal, 1, ImpulseResp ); 
   subplot( 2,2,3 );
   title('Convolution')
   plot( Conv );

   fftConv = fft(Conv);

   FIR = fft( ImpulseResp );

   Inverse = ifft( fftConv./FIR );

   subplot( 2,2,4 );
   plot( abs( Inverse ) );
  1. Как сказано, БПФ предполагает, что сигнал является периодическим. Если вы хотите анализировать апериодический сигнал, вам нужно использовать алгоритм «оконного БПФ» для улучшения качества. Это почти то же самое, но вы должны умножить свой ввод на специальную оконную функцию (например, Блэкман). Добавление отступов также работает, но это наиболее ресурсоемкий способ.
person Dmitry Galchinsky    schedule 11.03.2013
comment
+1: Учитывая тот факт, что вопрос касается conv, вместо filter я бы предложил сохранить conv и дополнить результаты БПФ нулями, используя команды fft(Conv, N) и fft(ImpulseResp, N), где N = length(TrueSignal) + length(ImpulseResp) - 1. Также нет необходимости делить fft(Conv) и fft(ImpulseResp) на длину сингла. - person Eitan T; 12.03.2013
comment
Большое спасибо! Я попробую это, мне определенно нужно теоретическое чтение по БПФ. - person agronskiy; 12.03.2013