1D прямая и обратная свертка с помощью преобразования Фурье

[Adventor]

Постановка задачи:
Есть некоторый исходный сигнал $s(t)$ (дискретизованный и квантованный), есть некоторая функция отклика $r(t)$, нужно получить свертку $c(t) = s(t) \ast r(t) $, а после чего восстановить исходный сигнал при помощи обратной свертки $ds(t) = c(t) \ast^{-1} r(t)$. Для решения задачи предлагается использовать БПФ.

Графическая иллюстрация ( взято из https://www.originlab.com )


Вопросы

1. Если сигнал имеет количество отсчетов $N$, то и БПФ также будет иметь $N$ отсчетов. В реализациях БПФ иногда встречается ограничение что $N$ должно быть кратно степени двойки. Я правильно понимаю, что в этом случае массив, к которому будет применено ПБФ, нужно дополнить нулями в конце так чтобы его длина стала кратна степени 2?

2. Таким же образом следует дополнить массив функции отклика $r(t)$, но уже до длины сигнала, увеличенной до кратности 2ке?

3. Чтобы получить свертку, нужно от поэлементного произведения массивов БПФ взять обратное БПФ, и потом от него взять действительную часть?

[Adventor]

Пока получается не очень. https://i.imgur.com/x3mOWsC.png
1. Свертка почему-то сдвинута по времени.
2. При вариации стандартного отклонения функции отклика сглаживания почему-то не происходит.
3. При вычислении обратной свертки присутствуют малые (околонулевые) значения, деление со всеми вытекающими.

Вот еще пример https://i.imgur.com/Gii0U9z.png
Есть шанс что прямая свертка считается правильно? Период на свертке растянулся в два раза.

[profrotter]

Adventor в сообщении #1298960 писал(а):

Пока получается не очень

А как должно быть очень? Что Вы вообще хотите получить?

[Adventor]

profrotter в сообщении #1299169 писал(а):

Adventor в сообщении #1298960 писал(а):

Пока получается не очень

А как должно быть очень? Что Вы вообще хотите получить?

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

[Adventor]

Кое-что в итоге получилось. Но возникли вопросы.

Случай первый: прямоугольный импульс и гаусс. Всё получилось. https://i.imgur.com/7S4SOg9.png

Случай второй: синус и гаусс. Свертка странная (но посчитана она вероятно правильно) и в качестве восстановленного сигнала шум. https://i.imgur.com/WS3xc59.png

Вопрос:
Почему во втором случае такая ситуация? Моё предположение, что я мог нарушить условия применимости. Но что именно я не понимаю.

Спасибо.

UPD:
Уменьшую длину сигнала, всё восстанавливается. Значит, я действительно что-то нарушил.
Либо уменьшать длину сигнала, либо увеличивать длину Фурье.

[arseniiv]

Если идти с самого начала, тут главное не забывать, что дискретное преобразование Фурье — это преобразование функции на «дискретной окружности», естественно такие функции считать периодическими, и свёртка с использованием ДПФ — это т. н. циклическая свёртка. Потому если нужна обычная свёртка двух непериодических функций, без достаточного количества нулей при свёртке левый кусок результата наложится на правый, и это будет уже не то. Другими словами, не стоит забывать, что свёртка последовательности $m$ значений и последовательности $n$ значений состоит из $m+n-1$ значений, так что надо место под это всё заранее выделить, беря размером БПФ не наименьшую степень двойки, не меньшую $\max\{m,n\}$, а наименьшую степень двойки, не меньшую $m+n-1$.

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

Если одна из функций периодическая, но период не укладывается в степень двойки, придётся делать более хитрые разрезания и сложения и, как и в первом случае (потому что «естественной» периодичностью мы воспользоваться не можем), отводить не меньше $m+n-1$ места под промежуточный результат.

[Pavia]

Решил тоже сделать.



-- Сб мар 24, 2018 18:21:13 --

Adventor в сообщении #1299185 писал(а):

Уменьшую длину сигнала, всё восстанавливается. Значит, я действительно что-то нарушил.
Либо уменьшать длину сигнала, либо увеличивать длину Фурье.

Я уменьшаю и увеличиваю у меня нормально восстанавливается.

-- Сб мар 24, 2018 18:33:08 --

Adventor

Adventor в сообщении #1298918 писал(а):

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

Скорее нет чем да. Есть несколько подходов зеркально отразить циклически продолжить экстраполировать сигнал на конце.
Думаю если взял среднее значение и заполнил им, то будет лучше чем если просто 0.

Adventor в сообщении #1298918 писал(а):

Чтобы получить свертку, нужно от поэлементного произведения массивов БПФ взять обратное БПФ, и потом от него взять действительную часть?

Да так не забываем что элементы Фурье комплексные числа и произведение соответственно должно быть комплексным.

Adventor в сообщении #1298960 писал(а):

1. Свертка почему-то сдвинута по времени.
Есть шанс что прямая свертка считается правильно? Период на свертке растянулся в два раза.

Это у вас линейная свёртка со сдвигом, а Фурье она циклическая.

Так что при таком подходе у вас будут не стыковки.

[Adventor]

arseniiv в сообщении #1299251 писал(а):

Если идти с самого начала, тут главное не забывать, что дискретное преобразование Фурье — это преобразование функции на «дискретной окружности», естественно такие функции считать периодическими, и свёртка с использованием ДПФ — это т. н. циклическая свёртка. Потому если нужна обычная свёртка двух непериодических функций, без достаточного количества нулей при свёртке левый кусок результата наложится на правый, и это будет уже не то. Другими словами, не стоит забывать, что свёртка последовательности $m$ значений и последовательности $n$ значений состоит из $m+n-1$ значений, так что надо место под это всё заранее выделить, беря размером БПФ не наименьшую степень двойки, не меньшую $\max\{m,n\}$, а наименьшую степень двойки, не меньшую $m+n-1$.

Спасибо за замечание. Длина Фурье получается действительно должна быть больше длины свертки $m+n-1$.

-- Вс мар 25, 2018 15:51:08 --

Pavia в сообщении #1299462 писал(а):

Это у вас линейная свёртка со сдвигом, а Фурье она циклическая.

Так что при таком подходе у вас будут не стыковки.

Спасибо за ответ. Цикличность свертки я в своем втором посте в этой теме уже учел.

-- Вс мар 25, 2018 15:53:09 --

Здесь же я могу рискнуть задать стыдный вопрос:
Мне казалось что Фурье от синусоиды будет содержать одну единственную частотную компоненту. На самом деле это очевидно не так. Почему? Это связано с дискретизацией сигнала?

[Pavia]

Adventor в сообщении #1299644 писал(а):

Здесь же я могу рискнуть задать стыдный вопрос:
Мне казалось что Фурье от синусоиды будет содержать одну единственную частотную компоненту. На самом деле это очевидно не так. Почему? Это связано с дискретизацией сигнала?

Глава 5 Растекание спектра в книге А.Б.Сергиенко Цифровая обработка сигналов.
Согласно теореме ряда Фурье нужен периодический сигнал с периодом $T$. Если частота синусоиды не совпадает, то мы отрезаем от неё кусок. И получается, что сигнал становится непериодическим. Это хорошо заметно, если соединить начальный и конечный кусок.

[Adventor]

Pavia в сообщении #1299682 писал(а):

Adventor в сообщении #1299644 писал(а):

Здесь же я могу рискнуть задать стыдный вопрос:
Мне казалось что Фурье от синусоиды будет содержать одну единственную частотную компоненту. На самом деле это очевидно не так. Почему? Это связано с дискретизацией сигнала?

Глава 5 Растекание спектра в книге А.Б.Сергиенко Цифровая обработка сигналов.
Согласно теореме ряда Фурье нужен периодический сигнал с периодом $T$. Если частота синусоиды не совпадает, то мы отрезаем от неё кусок. И получается, что сигнал становится непериодическим. Это хорошо заметно, если соединить начальный и конечный кусок.

К сожалению, не заезжаю...

Пример.
1 секунда времени и синусоида частотой 1Гц.
Частота сэмплирования 150Гц.
Отсчетов в сигнале 151.
Длина Фурье 1024 отсчета.

Итого один период периодического сигнала. Фурье спектр не простой.
https://i.imgur.com/evDe8Xv.png

Если сделать длину одного периода синусоиды равной длине Фурье, в частотном спектре получится один ненулевой отсчет.

UPD:
Синусоида дополняется нулями в конце чтобы её длина была равна длине Фурье. Значит, сигнал получается уже не синусоида. И поэтому Фурье спектр будет содержать уже не одну компоненту. Правильно?

Но периодичность сигнала же сохраняется?

Получается, в данном случае ради ДПФ мы меняем исходный сигнал. И получаем ДПФ уже другого, измененного сигнала? Это вообще правомерно?

[profrotter]

Adventor в сообщении #1299717 писал(а):

Получается, в данном случае ради ДПФ мы меняем исходный сигнал. И получаем ДПФ уже другого, измененного сигнала? Это вообще правомерно?

Не ради ДПФ, а ради БПФ. Но чаще вообще ради того, чтобы иметь тот шаг дискретизации спектра, который мы хотим иметь. Дополнение нулями исходного дискретного сигнала никак не влияет на его спектр, который определяется по известной формуле $S_d(\omega)=\sum_{n=0}^{N-1}s[n]e^{-j\omega nT}$. Если дописали к сигналу из $N$ отсчётов ещё $M$ нулевых отсчётов, то по факту это никак не изменяет верхний предел записанной суммы, поскольку суммировать нулевые слагаемые можно сколько угодно. Однако ДПФ в обоих случаях будет отличаться объёмом. В первом случае (до дополнения нулями) это будет ДПФ объёмом $N$, а во втором случае - ДПФ объёмом $N+M$. Но ДПФ, как известно, даёт отсчёты спектра дискретного сигнала $S_d(\omega)$, взятые с шагом $\Omega =\frac{2\pi}{NT}$. Соответственно при объёме ДПФ $N+M$ шаг дискретизации спектра будет $\Omega_{\text{с нулями}} =\frac{2\pi}{(N+M)T}<\Omega$.

Что касается одного периода синуса, то один период его спектра отнюдь не одна спектральная линия. Но при ДПФ спектр оказывается дискретизирован так, что только один отсчёт спектра отличается от нуля, а все остальные отсчёты приходятся в нулевые значения спектра. Если добавляются нулевые отсчёты, то спектр дискретизируется гуще и становятся видны отличные от нуля отсчёты спектра.

Этот эффект проще наблюдать даже не на гармонике, а на прямоугольном импульсе. Если применить ДПФ к $N$ единицам, то получим одну спектральную линию в нуле, а остальные значения ДПФ на периоде нулевыми. Если добавить нулевых отсчётов, то увидим дискретный спектр, похожий на отсчёты спектра прямоугольного импульса.

При желании вот здесь есть программа, которая показывает все спектры в накладку http://circuits-signals.narod.ru/DSpectr.zip

[Pavia]

Этот парадокс выше моего понимания. Дальше я пишу свои мысли они не закончены.

Если мы из спектра выберем определенные частоты, то нам это позволит создать однозначно определённую систему из коэффициентов.
И в ней будет строго одна частота.
Про то что там(не в спектре а в системе) одна частота описано у Фихтенгольц Г.М.-Основы математического анализа. Том 2-Наука (1968) пункт 397 страницы 374-376
Если нужно это свойство, то добавлять нули нет смысла.

Adventor в сообщении #1299717 писал(а):

Получается, в данном случае ради ДПФ мы меняем исходный сигнал. И получаем ДПФ уже другого, измененного сигнала?

Да это другой сигнал он перестал быть ортогональным к синусам и косинусам у него другой период он перестал быть симметричным. Или другое объяснение если продлить эту синусоиду, то концы её не сойдутся.

Можем мы применять? Несмотря на то, что все свойства нарушены, да можем применять. Во-первых остались какие-то свойство не нарушены к примеру интегрируемость. Во-вторых это же не усечении "информации", а расширение. Вот эту лишнюю прибавку мы всегда можем убрать и восстановить исходный сигнал. Нарушенные свойства очень сильно мешают. Поэтому инженеры и математики всячески стараются компенсировать и устранить их. Тежи нули можно добавить с двух сторон свёртку делят на казуальную и антиказульную часть и другие методы. Как то отказаться от свёртке в Фурье области, а реализовать её прямыми методами хотя там свои подводные камни.

[Adventor]

profrotter в сообщении #1299737 писал(а):

Дополнение нулями исходного дискретного сигнала никак не влияет на его спектр, который определяется по известной формуле $S_d(\omega)=\sum_{n=0}^{N-1}s[n]e^{-j\omega nT}$. Если дописали к сигналу из $N$ отсчётов ещё $M$ нулевых отсчётов, то по факту это никак не изменяет верхний предел записанной суммы, поскольку суммировать нулевые слагаемые можно сколько угодно. Однако ДПФ в обоих случаях будет отличаться объёмом. В первом случае (до дополнения нулями) это будет ДПФ объёмом $N$, а во втором случае - ДПФ объёмом $N+M$. Но ДПФ, как известно, даёт отсчёты спектра дискретного сигнала $S_d(\omega)$, взятые с шагом $\Omega =\frac{2\pi}{NT}$. Соответственно при объёме ДПФ $N+M$ шаг дискретизации спектра будет $\Omega_{\text{с нулями}} =\frac{2\pi}{(N+M)T}<\Omega$.

Не понимаю, вот здесь у меня не сходится.

Случай первый.
Синусоида 4 Герца. 1 секунда. Число отсчетов в синусоиде 1024. Длина сигнала 1024 отсчетов. Длина Фурье 1024 отсчетов.

Результат -- единственный ненулевой отсчет в спектре четко выделяется единственная гармоника.
https://i.imgur.com/7KqsGHo.png

Случай второй
Та же самая синусоида 4 Герца. 1 секунда. Число отсчетов в синусоиде 1024. Но длина сигнала 2048 отсчетов. Длина Фурье 2048 отсчетов.
Т.е. к синусоиде добавлены нули.

Результат -- спектр из многих гармоник.
https://i.imgur.com/Qw1lHk1.png

Так что может значить "дополнение нулями исходного дискретного сигнала никак не влияет на его спектр"? Если часть сигнала представлять в виде нулей, это будет и совершенно другой сигнал и соответственно другой спектр. Разве не так?

[profrotter]

Adventor в сообщении #1300835 писал(а):

Так что может значить

Это может значить то, что предложения следует дочитывать до конца:

profrotter в сообщении #1299737 писал(а):

Дополнение нулями исходного дискретного сигнала никак не влияет на его спектр, который определяется по известной формуле $S_d(\omega)=\sum_{n=0}^{N-1}s[n]e^{-j\omega nT}$.

В приведённую формулу можете подставлять сколько угодно новых нулевых значений сигнала и её результат никак не изменится, поскольку прибавление нуля не изменяет сумму. Это очевидно и не требует каких-либо иллюстраций. Здесь сам дискретный сигнал из $N$ отсчётов рассматривается как непериодический, длительностью $T_c=(N-1)T\approx NT$. (Единичкой при больших $N$ пренебрегают.) Период дискретизации спектра в этом случае можно выбрать в соответствии с теоремой Котельникова в частотной области как $\Omega\leq \frac{2\pi}{T_c}=\frac{2\pi}{NT} $. Возьмём максимальный период дискретизации спектра $\Omega = \frac{2\pi}{NT}$ и получим его отсчёты в виде: $S_d(k\Omega)=\sum_{n=0}^{N-1}s[n]e^{-j\frac{2\pi}{N}nk}$. ДПФ, таким образом, даёты отсчёты спектра дискретного сигнала, взятые с максимальным шагом.

Добавил 19:54 мск 01.04.2018

Предполгаю, что полезно будет просмотреть вот этот материал http://circuits-signals.narod.ru/part25.pdf (по ссылке pdf - лекция "Спектральный анализ дискретных сигналов"). Всё о чём там сказано я кратко описал, но там подробнее.
Пункты 5.4 и 5.6 особенно.

[Pavia]

Adventor

Adventor в сообщении #1300835 писал(а):

Так что может значить "дополнение нулями исходного дискретного сигнала никак не влияет на его спектр"? Если часть сигнала представлять в виде нулей, это будет и совершенно другой сигнал и соответственно другой спектр. Разве не так?

Это оптический обман зрения. Дело в том что у вас спектр нарисован непрерывной линией когда, как для дискретного сигнала он будет дискретным. В виде точек иногда рисуют в в виде стрелочек. А вот между стрелочками он может быть каким угодно. Когда вы поменяли размер преобразования Фурье вы увидели, то что до этого не видели- появились промежуточные точки.
Иногда это называют интерполяций сигнала. Как известно прямую через 2 точки мы можем провести одним способом, а кривую сколь угодно бесконечным числом способов. Вот это я и называю обманом зрения.

Осталось понять правильная эта интерполяция или нет. И какая правильная? Что это интерполяция приводит к ошибкам и так очевидна. Хотя profrotter и продолжает упрямствовать.

Теорема Котельникова говорит что она может восстановить гармонический сигнал. Она не применима не к гармоническим сигналам. Гармоничность у вас нарушена и это очевидно и применять вы её не в праве. Более того лучше читать Найквиста, который говорит, что для неорганических сигналов спектр должен идти до бесконечности.

Но любой не гармонический сигнал можно приближёно описать гармоническим. Да можно.

Но вот для дискретного сигнала у нас уже потерна информация между его отсчётами. Я говорю про дискретность сигнала. Для того чтобы правильно сделать интерполяцию надо априори знать информацию о сигнале. - profrotter пренебрегает этим.

В большинстве случаев она недоступна. Поэтому для качественной оценки мы можем дополнить сигнал нулями. А можем и просто линейно интерполировать. Разница невелика.

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