Почему результаты свертки имеют разную длину при выполнении во временной области vs в частотной области?

Я не эксперт DSP, но я понимаю, что есть два способа, которыми я могу применить дискретный фильтр во временной области к форме дискретной временной области. Во-первых, они свертывают их во временной области, а во-вторых, принимают БПФ обоих, умножают оба сложных спектра и принимают IFFT результата. Одним из ключевых различий в этих методах является второй подход, подлежащий циклической свертке.

В качестве примера, если фильтр и формы осциллограмм имеют длительность N точек, первый подход (т.е. свертка) дает результат, который имеет длину N + N-1, где первая половина этого ответа заполняется фильтром а вторая половина - опорожнение фильтра. Чтобы получить устойчивый ответ, фильтр должен иметь меньше точек, чем фильтруемый сигнал.

Продолжая этот пример со вторым подходом и предполагая, что данные формы дискретной временной области все являются реальными (не сложными), БПФ фильтра и формы волны создают БПФ из N точек в длину. Умножение обоих спектров. IFFT-результат приводит к тому, что результат во временной области также содержит N точек. Здесь ответ, когда фильтр заполняется и опустошается, перекрывает друг друга во временной области, и нет ответа устойчивого состояния. Это эффект круговой свертки. Чтобы этого избежать, обычно размер фильтра будет меньше размера осциллограммы, и оба будут иметь нулевую ширину, чтобы пространство для свертки частоты расширялось во времени после IFFT произведения двух спектров.

Мой вопрос заключается в том, что я часто вижу работу в литературе у хорошо зарекомендовавших себя специалистов/компаний, где они имеют дискретную (реальную) формулу во временной области (N точек), они FFT ее, умножают на некоторый фильтр (также N точек) и IFFT результат для последующей обработки. Мое наивное мышление заключается в том, что этот результат не должен содержать устойчивого ответа и, следовательно, должен содержать артефакты из заполнения/опорожнения фильтра, что приведет к ошибкам в интерпретации полученных данных, но я должен что-то упустить. При каких обстоятельствах это может быть обоснованным подходом?

Любое понимание было бы оценено

Ответ 1

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

Вернуться к первому абзацу: Обычно в DSP самая большая проблема, с которой я сталкиваюсь с БПФ, заключается в том, что они не являются причинами, и по этой причине я часто предпочитаю оставаться во временной области, используя, например, FIR и IIR-фильтры.

Update:

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

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

FFT. В подходе БПФ/свертки можно было бы получить образец, скажем, 256 точек, и сверлить его с помощью wfm, который является постоянным для первой половины и будет равен нулю во втором половина. Вопрос здесь (даже после того, как эта система провела несколько циклов), что определяет значение первой точки результата? FFT предполагает, что wfm является круговым (т.е. бесконечно-периодическим), так: либо первая точка результата определяется последними 127 (т.е. будущими) образцами wfm (пропуская середину wfm), либо 127 нулями если вы ноль-pad. Неправильно.

FIR. Другой подход - реализовать среднее значение с помощью FIR-фильтра. Например, здесь можно было бы использовать среднее значение значений в очереди регистра FIFO 128. То есть, по мере того, как каждый образец точки приходит, 1) помещает его в очередь, 2) удаляет старейший элемент, 3) усредняет все 128 элементов, оставшихся в очереди; и это ваш результат для этого примера. Этот подход работает непрерывно, обрабатывая по одной точке за раз и возвращая отфильтрованный результат после каждого образца и не имеет ни одной из проблем, которые возникают из БПФ, когда он применяется к конечным фрагментам выборки. Каждый результат - это только среднее значение текущего образца и 127 образцов, которые были перед ним.

В документе, который цитирует OP, применяется подход, гораздо более похожий на FIR-фильтр, чем на фильтр FFT (обратите внимание, что фильтр в документе более сложный, и вся бумага в основном представляет собой анализ этого фильтра.) См., например, эта бесплатная книга, в которой описывается, как анализировать и применять различные фильтры, а также отметить, что подход Лапласа к анализу FIR и IIR фильтры аналогичны тем, что найдено в цитируемой статье.

Ответ 2

Здесь приведен пример свертки без нулевого заполнения для ДПФ (круговая свертка) против линейной свертки. Это свертка последовательности M = 32 с длиной L = 128 (с использованием Numpy/Matplotlib):

f = rand(32); g = rand(128)
h1 = convolve(f, g)
h2 = real(ifft(fft(f, 128)*fft(g)))
plot(h1); plot(h2,'r')
grid()

Первые точки M-1 отличаются друг от друга, и они сокращаются по точкам M-1, поскольку они не были заполнены нулями. Эти различия являются проблемой, если вы выполняете блок-свертку, но для преодоления этой проблемы используются такие методы, как перекрытие и сохранение или перекрытие и добавление. В противном случае, если вы просто вычисляете одноразовую операцию фильтрации, действительный результат будет начинаться с индекса M-1 и заканчиваться индексом L-1 с длиной L-M + 1.

Что касается цитируемой статьи, я просмотрел их код MATLAB в приложении A. Я думаю, что они допустили ошибку в применении функции передачи Hfinal к отрицательным частотам без ее первого сопряжения. В противном случае на их графиках видно, что дрожание часов является периодическим сигналом, поэтому использование круговой свертки отлично подходит для стационарного анализа.

Edit: Что касается сопряжения передаточной функции, PLL имеют действительный импульсный отклик, и каждый вещественнозначный сигнал имеет сопряженный симметричный спектр. В коде вы можете видеть, что они просто используют Hfinal [N-i] для получения отрицательных частот без принятия сопряженного. Я построил их функцию передачи от -50 МГц до 50 МГц:

N = 150000                    # number of samples. Need >50k to get a good spectrum. 
res = 100e6/N                 # resolution of single freq point  
f = res * arange(-N/2, N/2)   # set the frequency sweep [-50MHz,50MHz), N points
s = 2j*pi*f                   # set the xfer function to complex radians 

f1 = 22e6       # define 3dB corner frequency for H1 
zeta1 = 0.54    # define peaking for H1 
f2 = 7e6        # define 3dB corner frequency for H2 
zeta2 = 0.54    # define peaking for H2    
f3 = 1.0e6      # define 3dB corner frequency for H3 

# w1 = natural frequency   
w1 = 2*pi*f1/((1 + 2*zeta1**2 + ((1 + 2*zeta1**2)**2 + 1)**0.5)**0.5)  
# H1 transfer function 
H1 = ((2*zeta1*w1*s + w1**2)/(s**2 + 2*zeta1*w1*s + w1**2))            

# w2 = natural frequency 
w2 = 2*pi*f2/((1 + 2*zeta2**2 + ((1 + 2*zeta2**2)**2 + 1)**0.5)**0.5)  
# H2 transfer function  
H2 = ((2*zeta2*w2*s + w2**2)/(s**2 + 2*zeta2*w2*s + w2**2))            

w3 = 2*pi*f3        # w3 = 3dB point for a single pole high pass function. 
H3 = s/(s+w3)       # the H3 xfer function is a high pass

Ht = 2*(H1-H2)*H3   # Final transfer based on the difference functions

subplot(311); plot(f, abs(Ht)); ylabel("abs")
subplot(312); plot(f, real(Ht)); ylabel("real")
subplot(313); plot(f, imag(Ht)); ylabel("imag")

Как вы можете видеть, реальный компонент имеет даже симметрию, а мнимый компонент имеет нечетную симметрию. В своем коде они вычисляли только положительные частоты для логарифмического графика (достаточно разумно). Однако для расчета обратного преобразования они использовали значения для положительных частот для отрицательных частот путем индексирования Hfinal [N-i], но забывали сопрягать его.

Ответ 3

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

Ответ 4

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

Как уже указывалось, FFT предполагает, что мы имеем бесконечный сигнал. Когда мы берем образец за конечное время T, это математически эквивалентно умножению сигнала с прямоугольной функцией.

Умножение во временной области становится сверткой в частотной области. Частотной характеристикой прямоугольника является синхронизирующая функция, то есть sin (x)/x. X в числителе является кикером, поскольку он умирает O (1/N).

Если у вас есть частотные составляющие, которые точно кратно 1/T, это не имеет значения, так как функция синхронизации равна нулю во всех точках, кроме той частоты, где она равна 1.

Однако, если у вас есть синус, который падает между 2 точками, вы увидите функцию синхронизации, отобранную на частотной точке. Он похож на увеличенную версию функции синхронизации, а сигналы "призрак", вызванные сверткой, сжимаются с 1/N или 6 дБ/октавой. Если у вас есть сигнал 60 дБ выше уровня шума, вы не увидите шума на 1000 частот слева и справа от основного сигнала, он будет зависеть от "юбок" функции синхронизации.

Если вы используете другое временное окно, вы получаете разную частотную характеристику, например, косинус умирает с 1/x ^ 2, есть специальные окна для разных измерений. Окно Ханнинга часто используется как окно общего назначения.

Дело в том, что прямоугольное окно, используемое при отсутствии какой-либо "функции окна", создает гораздо худшие артефакты, чем хорошо выбранное окно. то есть путем "искажения" временных выборок мы получаем гораздо лучшую картину в частотной области, которая ближе напоминает "реальность", или, скорее, "реальность", которую мы ожидаем и хотим видеть.