[БендеровецЪ]

Подумалось что может быть полезным.
Как сделать коректирующий FIR из данных ипортированых, например, из ЛТ спайс (ну или откуда-то еще) через ifft. Как всегда можно воспользоваться специальными програмками для этого. Но иногда задача немного нестандартна по какой-либо причине, что делает их неудобными. И вобще меня бесит когда я не вижу что эти програмки делают внутри :)

Вобщем я делаю это в numpy где-то так:

Импортировать необходимые библиотеки:

Код:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy import signal
from scipy.io import wavfile
import time as tm
% matplotlib inline

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

Код:

data_raw = np.loadtxt("example__complx.txt", delimiter=",")

Задать fs:

Код:

fs = 96e3

Разделить колонки на частоты и коплексное число коэф. передачи:

Код:

freq = data_raw[:,0]
resp = data_raw[:,1] + 1j*data_raw[:,2]

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

Задать длинну фильтра и получить сетку частот (вобщем-то не должно быть 2^n, просто тогда это будте dft а не fft)

Код:

fft_length = np.power(2, 12)
fft_freq = np.fft.rfftfreq(fft_length*2)*fs

Проинтерполировать АФЧХ в новую сетку частот:

Код:

ifft_input = (np.interp(fft_freq, freq, data_raw[:,1])) + 1j*(np.interp(fft_freq, freq, data_raw[:,2]))

Где-то на этом этапе можно добавить сглаживание, если необходимо. nympy знает всяких там савитцких-голе лично, соотв смотрет одноименную функцию.

Поскольку меня интересует тлько корекция ФЧХ, то нормализирую комплесные числа до плоской АЧХ (вобще можно записать сильно по разному, например просто поделить на абсолютное значение):

Код:

ifft_input_flat = np.cos(np.angle(ifft_input)) - 1j*np.sin(np.angle(ifft_input))

минус перед 1j потому что нас интересует обратная фаза. Если нужна кореция АЧХ, то можно просто взять ресипрокал по амплитуде.

выполнить ifft. Поскольку последовательность должна быть реальной, то можно воспользоваться пункцией irfft, которая подразумевает хермитиан симметрию входной последовательности, поэтому можно обойтись только "половиной" входных данных:

Код:

fir = np.roll(np.fft.irfft(ifft_input_flat, n = fft_length*2), fft_length)

В этой же строчке я делаю roll ("прокрутку") на пол длинны окна, т.к. fft "центруется" относительно начала подразумевая циклический сигнал.

Нормализировать фильтр к Ку=1

fir /= np.sum(fir)

Контроль того что получилось на картинке:

Код:

plt.figure()
plt.plot(freq, np.unwrap(np.angle(resp)))
plt.plot(fft_freq, np.unwrap(np.angle(ifft_input_flat)), )
plt.xscale("log")
plt.show()

plt.figure()
plt.plot(fir)
plt.show()

Сохранить все это в wav файл (ну или в каком-там надо формате):

Код:

sp.io.wavfile.write('filename.wav', fs, np.float32(fir))

Готово. Все просто и прозрачно.

Для контроля можно вычитать данные из файла

Код:

fs, sig = wavfile.read('filename.wav')

И отобразить на графике все это дело:

Код:

plt.figure()
plt.plot((sig))
plt.show()