Как реализовать фильтр, например scipy.signal.lfilter

Question

Я сделал прототип на python, который я конвертирую в приложение iOS. К сожалению, все приятные черты scipy и numpy недоступны в объективе C. Поэтому, по-видимому, мне нужно реализовать фильтр в объективе C с нуля. В качестве первого шага я пытаюсь реализовать IIR с нуля в python. Если я смогу понять, как это сделать в python, я смогу закодировать его в C.

В качестве примечания, я был бы признателен за любые предложения по ресурсам по фильтрации в iOS. Как новичок в объективе C, который привык к matlab и python, я шокирован тем, что такие вещи, как Audio Toolboxes и Accelerate Frameworks и Amazing Audio Engines, не имеют эквивалента scipy.signal.filtfilt, а также функций дизайна фильтров, таких как scipy. signal.butter и т. д.

Итак, в следующем коде я реализую фильтр пятью способами. 1) scipy.signal.lfilter (для сравнения), 2) форму пространства состояний с использованием матриц A, B, C, D, рассчитанных по масляной функции Matlab. 3) форму пространства состояний с использованием матриц A, B, C, D, рассчитанных scipy.signal.tf2ss. 4) Прямая форма I, 5) Прямая форма II.

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

import numpy as np 
from scipy.signal import butter, lfilter, tf2ss 

# building the test signal, a sum of two sines; 
N = 32 
x = np.sin(np.arange(N)/6. * 2 * np.pi)+\ 
    np.sin(np.arange(N)/32. * 2 * np.pi) 
x = np.append([0 for i in range(6)], x) 

# getting filter coefficients from scipy 
b,a = butter(N=6, Wn=0.5) 

# getting matrices for the state-space form of the filter from scipy. 
A_spy, B_spy, C_spy, D_spy = tf2ss(b,a) 

# matrices for the state-space form as generated by matlab (different to scipy's) 
A_mlb = np.array([[-0.4913, -0.5087, 0, 0, 0, 0], 
     [0.5087, 0.4913, 0, 0, 0, 0], 
     [0.1490, 0.4368, -0.4142, -0.5858, 0, 0], 
     [0.1490, 0.4368, 0.5858, 0.4142, 0, 0], 
     [0.0592, 0.1735, 0.2327, 0.5617, -0.2056, -0.7944], 
     [0.0592, 0.1735, 0.2327, 0.5617, 0.7944, 0.2056]]) 

B_mlb = np.array([0.7194, 0.7194, 0.2107, 0.2107, 0.0837, 0.0837]) 

C_mlb = np.array([0.0209, 0.0613, 0.0823, 0.1986, 0.2809, 0.4262]) 

D_mlb = 0.0296 

# getting results of scipy lfilter to test my implementation against 
y_lfilter = lfilter(b, a, x) 

# initializing y_df1, the result of the Direct Form I method. 
y_df1 = np.zeros(6) 

# initializing y_df2, the result of the Direct Form II method. 
# g is an array also used in the calculation of Direct Form II 
y_df2 = np.array([]) 
g = np.zeros(6) 

# initializing z and y for the state space version with scipy matrices. 
z_ss_spy = np.zeros(6) 
y_ss_spy = np.array([]) 

# initializing z and y for the state space version with matlab matrices. 
z_ss_mlb = np.zeros(6) 
y_ss_mlb = np.array([]) 


# applying the IIR filter, in it's different implementations 
for n in range(N): 
    # The Direct Form I 
    y_df1 = np.append(y_df1, y_df1[-6:].dot(a[:0:-1]) + x[n:n+7].dot(b[::-1])) 

    # The Direct Form II 
    g = np.append(g, x[n] + g[-6:].dot(a[:0:-1])) 
    y_df2 = np.append(y_df2, g[-7:].dot(b[::-1])) 

    # State space with scipy's matrices 
    y_ss_spy = np.append(y_ss_spy, C_spy.dot(z_ss_spy) + D_spy * x[n+6]) 
    z_ss_spy = A_spy.dot(z_ss_spy) + B_spy * x[n+6] 

    # State space with matlab's matrices 
    y_ss_mlb = np.append(y_ss_mlb, C_mlb.dot(z_ss_mlb) + D_mlb * x[n+6]) 
    z_ss_mlb = A_mlb.dot(z_ss_mlb) + B_mlb * x[n+6] 


# getting rid of the zero padding in the results 
y_lfilter = y_lfilter[6:] 
y_df1 = y_df1[6:] 
y_df2 = y_df2[6:] 


# printing the results 
print "{}\t{}\t{}\t{}\t{}".format('lfilter','matlab ss', 'scipy ss', 'Direct Form I', 'Direct Form II') 
for n in range(N-6): 
    print "{}\t{}\t{}\t{}\t{}".format(y_lfilter[n], y_ss_mlb[n], y_ss_spy[n], y_df1[n], y_df2[n]) 

И выход:

lfilter   matlab ss  scipy ss  Direct Form I Direct Form II 
0.0    0.0    0.0    0.0    0.0 
0.0313965294015 0.0314090254837 0.0313965294015 0.0313965294015 0.0313965294015 
0.225326252712 0.22531468279 0.0313965294015 0.225326252712 0.225326252712 
0.684651781013 0.684650012268 0.0313965294015 0.733485689277 0.733485689277 
1.10082931381 1.10080090424 0.0313965294015 1.45129994748 1.45129994748 
0.891192957678 0.891058879496 0.0313965294015 2.00124367289 2.00124367289 
0.140178897557 0.139981099035 0.0313965294015 2.17642377522 2.17642377522 
-0.162384434762 -0.162488434882 0.225326252712 2.24911228252 2.24911228252 
0.60258601688 0.602631573263 0.225326252712 2.69643931422 2.69643931422 
1.72287292534 1.72291129518 0.225326252712 3.67851039998 3.67851039998 
2.00953056605 2.00937857026 0.225326252712 4.8441925268 4.8441925268 
1.20855478679 1.20823164284 0.225326252712 5.65255635018 5.65255635018 
0.172378732435 0.172080718929 0.225326252712 5.88329450124 5.88329450124 
-0.128647387408 -0.128763927074 0.684651781013 5.8276996139 5.8276996139 
0.47311062085 0.473146568232 0.684651781013 5.97105082682 5.97105082682 
1.25980235112 1.25982698592 0.684651781013 6.48492462347 6.48492462347 
1.32273336715 1.32261397627 0.684651781013 7.03788646586 7.03788646586 
0.428664985784 0.428426965442 0.684651781013 7.11454966484 7.11454966484 
-0.724128943343 -0.724322419906 0.684651781013 6.52441390718 6.52441390718 
-1.16886662032 -1.16886884238 1.10082931381 5.59188293911 5.59188293911 
-0.639469994539 -0.639296371149 1.10082931381 4.83744942709 4.83744942709 
0.153883055505 0.154067363252 1.10082931381 4.46863620556 4.46863620556 
0.24752293493 0.247568224184 1.10082931381 4.18930262192 4.18930262192 
-0.595875437915 -0.595952759718 1.10082931381 3.51735265599 3.51735265599 
-1.64776590859 -1.64780228552 1.10082931381 2.29229811755 2.29229811755 
-1.94352867959 -1.94338167159 0.891192957678 0.86412577159 0.86412577159 

Answer 1

Итак, я наконец-то нашел ту часть рамки ускорения я искал.

Я использовал фильтр в первую очередь для понижающей дискретизации; вам необходимо отфильтровать перед понижающей дискретизацией, чтобы избежать наложения псевдонимов. Accelerate's vDSP_zrdesamp - это функция, которую я хотел все время.

Кроме того, для одного только фильтрования ipodEQ аудиоустройство, как правило, правильный выбор: (подтип kAudioUnitSubType_AUiPodEQ)

Если вы на самом деле нужно реализовать фильтр нуля, форма пространства состояний кажется лучшим.

Все еще не ответило: почему бы мне не использовать мои формы I и II?

Answer 2

взгляд на FIR wiki, и эта формула:

поэтому сначала вы создаете окно Хэмминга самостоятельно (я все еще использую питона но вы можете перевести его на c/cpp):

def getHammingWin(n): 
    ham=[0.54-0.46*np.cos(2*np.pi*i/(n-1)) for i in range(n)] 
    ham=np.asanyarray(ham) 
    ham/=ham.sum() 
    return ham 

затем свой собственный фильтр нижних частот:

def myLpf(data, hamming): 

    N=len(hamming) 
    res=[] 
    for n, v in enumerate(data): 
     y=0 
     for i in range(N): 
      if n<i: 
       break 
      y+=hamming[i]*data[n-i] 
     res.append(y) 
    return np.asanyarray(res) 
    pass