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

Простенькие генераторчики используя numpy, варианты вывода в int16 и float32. Переделать int16 в int32 несложно редактированием одного параметра. Задание длинны и частоты тоже вполне очевидно.
Несложно переделать под любое количество частот/форму сигнала. Количество дизера тоже параметризировано.

ЗЫ. Форум так и не позволяет атачить фалы с расширением .ipynb

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

Вобще все эти дела с ффт во времени делается и в матлабе/пайтоне. Надо только приблизительно представлять что хочется сделать. Вот пример чирпа записаного марленом. Видно пачти все, и гармоники, и НЧ срань.

ЗЫ Помоему промахнулся с горизонтальной шкалой в два раза.

Понял почему - делал 50% перекрытие окна.

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

Подумалось что может быть полезным.
Как сделать коректирующий 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()

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

А вот пример как запилить коректирующий IIR по типу "линквица".

Для начала нада запилить фунцию конвертации аналогова фильтра в IIR. Там есть свои приколы, и даже несколько вариантов имплементации (импулс инвариант, билинейный).
Подробнее обсуждается сдесь: https://dsp.stackexchange.com/questions/...pole/19211
Примеры: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Ну и тут как всегда удобно сделать свою библиотечку (во вложении). Принимает три параметра: вектор b вектор a и 1/fs синтезируемого iir

Использование выглядит приблизительно вот так:

Код:

from analog2iir import analog_iir
b_z_cor, a_z_cor = analog_iir(a_s, a_s_new, 1./fs)

(для моего примера fs = 96000.)

Собственно предидущая строка уже содержит то как делается коректирующий фильтр для системы (b_s, a_s) для получения желаемого отклика аналогичному для системы с (b_s_new, a_s_new) прт одинаковых a.

Собсна если кто помнит как делается "коректор Линквица". В числитель ставится знаменатель имеющейся системы, а в знаменатель - знаменатель целевой. Что я и сделал сразу вписав "a_s" и "a_s_new" в соотв места.

Для проверки хорошо бы исследовать АЧХ и ФЧХ исходной системы, целефой системы и корректирующего фильтра. Удобно использоваю уже имеюшиеся функции freqs и freqz

Код:

w_s, h_s = sp.signal.freqs(b_s, a_s, worN=np.logspace(1, 5, 1000)*2.*np.pi)
w_s_new, h_s_new = sp.signal.freqs(b_s_new, a_s_new, worN=np.logspace(1, 5, 1000)*2.*np.pi)
w_z, h_z = signal.freqz(b_z_cor, a_z_cor, worN=np.logspace(1, 5, 1000)*2.*np.pi/fs)

строим все это безобразие на графике:

Код:

plt.figure()
plt.plot(w_s/(2*np.pi), np.absolute(h_s))
plt.plot(fs*w_z/(2*np.pi), np.absolute(h_z))
plt.plot(w_s/(2*np.pi), np.absolute(h_s_new))
plt.xscale('log')
plt.yscale('log')
plt.xlim(10, 1e5)
plt.ylim(1e-5, 1)
plt.grid(which = 'both')
plt.show()

plt.figure()
plt.plot(w_s/(2*np.pi), np.angle(h_s))
plt.plot(w_s/(2*np.pi), np.angle(h_s_new))
plt.plot(fs*w_z/(2*np.pi), np.angle(h_z))
plt.xscale('log')
plt.xlim(10, 1e5)
plt.grid(which = 'both')
plt.show()

и получаем



Для примера использовал следующие вводные:

Код:

fc = 500. #Hz
fc_new = 100. #Hz
w_s = 2.*np.pi*fc
w_s_new = 2.*np.pi*fc_new
q_s = 0.7
q_s_new = 0.7

b_s = np.array((np.square(1./w_s) , 0, 0))
a_s = np.array((np.square(1./w_s) , 1./(w_s*q_s) , 1))

b_s_new = np.array((np.square(1./w_s_new) , 0, 0))
a_s_new = np.array((np.square(1./w_s_new) , 1./(w_s_new*q_s_new) , 1))

При желании можно просимулировать сисемы во временной области (все эти ваши импульсные отклики или ступеньки) ипользуя готовые функции lsim и dlsim.

NB ну если где-то с чем-то промахнулся то просьба не пинать, это все-таки не моя специальность :)

[mellowman]

а тут есть IIR для спайса: http://ltspicegoodies.ltwiki.org/Filter_...mples.html для тех кто не любит питон :)

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

Для тех кто не любит пайтон есть октав :)