Jakie jest najlepsze przybliżenie IIR (filtr AR) pierwszego rzędu do filtra średniej ruchomej (filtr FIR)?


24

Załóżmy następujący filtr IIR pierwszego rzędu:

y[n]=αx[n]+(1α)y[n-1]

Jak mogę wybrać parametr IIR przybliża tak dobrze, jak to możliwe, FIR, który jest średnią arytmetyczną ostatnich próbek:αk

z[n]=1kx[n]+1kx[n-1]++1kx[n-k+1]

Gdzie , co oznacza, że ​​dane wejściowe dla IIR mogą być dłuższe niż a mimo to chciałbym uzyskać najlepsze przybliżenie średniej z ostatnich danych wejściowych.n[k,)kk

Wiem, że IIR ma nieskończoną odpowiedź impulsową, dlatego szukam najlepszego przybliżenia. Byłbym szczęśliwy dla rozwiązania analitycznego, czy to dla funkcji kosztu czy .L 1L2L1

Jak można rozwiązać problemy z optymalizacją, biorąc pod uwagę tylko IIR pierwszego rzędu.

Dzięki.


Czy musi dokładnie podążać ]? y[n]=αx[n]+(1α)y[n1]
Phonon

1
To z pewnością stanie się bardzo słabym przybliżeniem. Czy nie stać cię na coś więcej niż IIR pierwszego rzędu?
leftaroundabout

Możesz zredagować swoje pytanie, aby nie używać w znaczeniu dwóch różnych rzeczy, np. Drugie wyświetlane równanie może brzmieć , a możesz chcieć powiedzieć, jakie dokładnie jest twoje kryterium „jak najlepiej”, np. czy chcesz być możliwie małe dla wszystkich , lub być tak małe, jak to możliwe dla wszystkich . z [ n ] = 1y[n]| y[n]-z[n]| n| y[n]-z[n]| 2nz[n]=1kx[n]++1kx[nk+1]|y[n]z[n]|n|y[n]z[n]|2n
Dilip Sarwate

@Phonon, tak, to musi być IIR pierwszego rzędu. Kryteria są proste, wynik powinien być jak najbardziej zbliżony do średniej z ostatnich k danych wejściowych do systemu, gdzie n [ k , inf ] . Byłbym szczęśliwy widząc wynik w obu przypadkach. Chociaż zakładam, że rozwiązanie analityczne jest opłacalne tylko dla | y [ n ] - z [ n ] | 2 . y[n]kn[k,inf]|y[n]z[n]|2
Royi,

Odpowiedzi:


10

Nie ma rozwiązania analitycznego dla będącego skalarem (tak myślę). Oto skrypt, który daje a dla danego K . Jeśli potrzebujesz go online, możesz zbudować LUT. Skrypt znajduje rozwiązanie, które minimalizujeααK

0πdw|H1(jw)H2(jw)|2

gdzie to odpowiedź częstotliwościowa FIR, a H 2 to odpowiedź częstotliwościowa IIR.H1H2

Nie określiłeś żadnego zakresu dla K. Ale chcę tylko wyjaśnić, że następujący system jest równoważny z twoim średnim filtrem i ma tę samą złożoność obliczeniową i twoje pierwsze zamówienie IIR!

H(z)=1K1zK1z1

function a = find_a(K)

w = 0.0001:0.001:pi;
as = [-1:0.001:-0.001  0.001:0.001:1];

E = zeros(size(as));
for idx=1:length(as)
    fJ = J(w,as(idx),K);
    E(idx) = sum(fJ);
end

[Emin, indx] = min(E)
a = as(indx)

function f = J(w,a,K)
    num = 2*(2-a)*(1-cos(w*K)) + 2*(cos(w*(K-1)) - cos(w)) - 2*(1-a)*(cos(w)-cos(w*(K+1)));
    den = (2-a)^2 + 1 + (1-a)^2 + 2*(1-a)*cos(2*w) - 2*(2-a)^2*cos(w);
    f = -(a/K)*num./den;
    f = f+(1/K^2)*(1-cos(w*K))./(1-cos(w))+a^2./(1+(1-a)^2-2*(1-a)*cos(w));
end

end

@Drazick Jest względnie prosty. Dwa wyrażenia dla IIR i FIR są podłączone do całki. Kluczem do znalezienia alternatywnego wyrażenia dla filtra FIR jest rozpoznanie postępu / serii geometrycznej. Wszystkie szczegóły znajdziesz tutaj: en.wikipedia.org/wiki/Geometric_progression#Geometric_series . W skrypcie funkcja J oblicza wyrażenie pod znakiem integralnym.
niaren

@niaren Wiem, że to jest stary post, więc jeśli pamiętasz: jak wyprowadzono twoją funkcję „f”? Kodowałem podobną rzecz, ale używając złożonych funkcji transferu dla FIR (H1) i IIR (H2), a następnie wykonując sumę (abs (H1 - H2) ** 2). Porównałem to z twoją sumą (fj), ale otrzymałem różne wyniki. Myślałem, że zapytam, zanim przeorę matematykę.
Dom

@Dom That That dawno temu i naprawdę nie pamiętam. Chyba właśnie przeszedłem proces opracowywania . Nie pamiętam, jak zweryfikowałem wyrażenie. Nie mam nic przeciwko [H1(jω)H2(jω)][H1(jω)H2(jω)]
powtórzeniu

@niaren Cześć, próbowałem wyprowadzić twoje wyrażenie, ale utknąłem przy sumowaniu złożonych frakcji. Popełniłem błąd w moim kodzie ... wydaje się, że twoja funkcja daje wyniki proporcjonalne do sumy (abs (H1 - H2) ** 2), co wskazuje, że jest poprawna.
Dom

16

Miło dyskutuje się na ten temat w przetwarzaniu sygnałów osadzonych za pomocą architektury mikro sygnałów , mniej więcej między stronami 63 i 69 . Na stronie 63 zawiera wyprowadzenie dokładnego filtru rekurencyjnej średniej ruchomej (który niaren podał w swojej odpowiedzi ),

H(z)=1N1zN1z1.

Dla wygody w odniesieniu do poniższej dyskusji odpowiada to równaniu różnicy:

yn=yn1+1N(xnxnN).

Przybliżenie, które wprowadza filtr do określonej przez ciebie postaci, wymaga przyjęcia, że , ponieważ (i cytuję z str. 68 ) „ y n - 1 jest średnią z x n próbek”. To przybliżenie pozwala nam uprościć poprzednie równanie różnicy w następujący sposób:xnNyn1yn1xn

yn=yn1+1N(xnyn1)yn=yn11Nyn1+1Nxnyn=(11N)yn1+1Nxn.

Ustawienie , dochodzimy do twojej pierwotnej postaci,yn=αxn+(1-α)yn-1, co pokazuje, że pożądany współczynnik (w odniesieniu do tego przybliżenia) wynosi dokładnie1α=1Nyn=αxn+(1α)yn1 (gdzieNjest liczbą próbek). 1NN

Czy to przybliżenie jest w pewnym sensie „najlepsze”? Z pewnością jest elegancki. Oto jak porównuje się odpowiedź amplitudowa [przy 44,1 kHz] dla N = 3, a gdy N wzrasta do 10 (przybliżenie na niebiesko):

N = 3 N = [3,10]


Jak sugeruje odpowiedź Petera , zbliżenie filtra FIR za pomocą filtra rekurencyjnego może być problematyczne pod normą dotyczącą najmniejszych kwadratów. Obszerną dyskusję na temat ogólnego rozwiązania tego problemu można znaleźć w pracy dyplomowej JOS, Techniki projektowania filtrów cyfrowych i identyfikacji systemu za pomocą aplikacji na skrzypce . Opowiada się za użyciem Normy Hankela, ale w przypadkach, w których reakcja fazowa nie ma znaczenia, obejmuje on także Metodę Kopca, która może w tym przypadku działać dobrze (i stosuje normę ). Szeroki przegląd technik zawartych w pracy można znaleźć tutaj . Mogą dać inne interesujące przybliżenia.L2


Jest to „elegancki” sposób na powiedzenie czegoś o pamięci filtra IIR pierwszego rzędu. Jego pamięć odpowiada . Przyjrzę się innym metodom, o których wspomniałeś. Dzięki. 1α
Royi,

Czy możesz wyjaśnić intuicyjnie, dlaczego zgodnie z normą LS ( ) nie ma rozwiązania? L2
Royi,

Nie jestem pewien, czy istnieje rozwiązanie LS w tym przypadku, czy jeszcze nie, po prostu wiem, że zwykle występują problemy ze zbieżnością w przypadku ogólnego problemu „przybliżenia FIR na podstawie IIR”. Zaktualizuję w / więcej informacji, kiedy będę mieć szansę.
Datageist

Cóż, jeśli funkcja kosztu zaproponowana przez Petera (pierwsza) jest właściwa, istnieje rozwiązanie. Przynajmniej według moich obliczeń.
Royi,

Świetny. Jestem ciekawy, jak „heurystyczne” podejście wypada w porównaniu z czymś bardziej kanonicznym. 1/N
Datageist

16

OK, spróbujmy wyprowadzić to, co najlepsze: tak, że współczynnik zx[n-m]toα(1-α

y[n]=αx[n]+(1α)y[n1]=αx[n]+(1α)αx[n1]+(1α)2y[n2]=αx[n]+(1α)αx[n1]+(1α)2αx[n2]+(1α)3y[n3]
x[nm] .α(1α)m

Najlepsze przybliżenie średnich kwadratów zminimalizuje:

J(α)=m=0k1(α(1α)m1k)2+m=kα2(1α)2m=m=0k1(α2(1α)2m2kα(1α)m+1k2)+α2(1α)2km=0(1α)2m=α21(1α)2k1(1α)2+2αk1(1α)k1(1α)+α2(1α)2k1(1α)2+1k=α21(1α)2+2k(1(1α)k)+1k=α22αα2+2k(1(1α)k)+1k=α2α+2k(1(1α)k)+1k
because the FIR coefficients are zero for m>k1.

Next step is to take derivatives and equate to zero.


Looking at a plot of the derived J for K=1000 and α from 0 to 1, it looks like the problem (as I've set it up) is ill-posed, because the best answer is α=0.

enter image description here


I think there's a mistake here. The way it should be according to my calculations is:

J(α)=m=0k1(α(1α)m1k)2+m=kα2(1α)2m=m=0k1(α2(1α)2m2kα(1α)m+1k2)+α2(1α)2km=0(1α)2m=α21(1α)2k1(1α)22αk1(1α)k1(1α)+1k+α2(1α)2k1(1α)2

Simplifying it according to Mathematica yields:

J(α)=α2α+2(1α)k1k

Using the following code on MATLAB yields something equivalent though different:

syms a k;

expr1 = (a ^ 2) * ((1 - ((1 - a) ^ (2 * k))) / (1 - ((1 - a) ^ 2)));
expr2 = ((2 * a) / k) * ((1 - ((1 - a) ^ (k))) / (1 - (1 - a)));
expr3 = (1 / k);
expr4 = ((a ^ 2) * ((1 - a) ^ (2 * k))) / (1 - ((1 - a) ^ (2)));

simpExpr = simplify(expr1 - expr2 + expr3 + expr4);

J(α)=2α2k2(1α)k+1k

Anyhow, those functions do have minimum.


So let's assume that we really only care about the approximation over the support (length) of the FIR filter. In that case, the optimization problem is just:

J2(α)=m=0k1(α(1α)m1k)2

Plotting J2(α) for various values of K versus α results in the date in the plots and table below.

For K = 8. αmin = 0.1533333
For K = 16. αmin = 0.08
For K = 24. αmin = 0.0533333
For K = 32. αmin = 0.04
For K = 40. αmin = 0.0333333
For K = 48. αmin = 0.0266667
For K = 56. αmin = 0.0233333
For K = 64. αmin = 0.02
For K = 72. αmin = 0.0166667

enter image description here

The red dashed lines are 1/K and the green lines are αmin, the value of α that minimizes J2(α) (chosen from alpha=[0:.01:1]/3;).


1
Was just going to post the exact same thing = )
Phonon

@Phonon: Feel free to continue! I've marked it as community wiki for that purpose.
Peter K.

The derivative w.r.t α is a series with an infinite number of terms (i.e. not a polynomial) that you have to set equal to 0 and then solve for α, and so some care (or possibly approximation) is going to be necessary.
Dilip Sarwate

Can someone please check and/or correct my working? :-)
Peter K.

@DilipSarwate, What would be the best approximation? Thanks.
Royi


3

I stumbled upon this old question and I would like to share my solution. As mentioned in other answers, there is no analytical solution, but the function to be minimized behaves nicely and the optimal value of α can be found easily with a few Newton iterations. There is also a formula to check the optimality of the result.

The impulse response of the length N FIR moving average filter is given by

(1)hFIR[n]=1N(u[n]u[nN])

where u[n] is the unit step function. The first order IIR filter

(2)y[n]=αx[n]+(1α)y[n1]

has the impulse response

(3)hIIR[n]=α(1α)nu[n]

The goal is now to minimize the squared error

(4)ϵ=n=0(hFIR[n]hIIR[n])2

Using (1) and (3), the error can be written as

ϵ(α)=n=0N1(α(1α)n1N)2+n=Nα2(1α)2n=α2n=0(1α)2n2αNn=0N1(1α)n+n=0N11N2=α21(1α)22αN1(1α)N1(1α)+1N(5)=α2α2N(1(1α)N)+1N,0<α<2

This expression is very similar to the one given in this answer, but it's not identical. The restriction on α in (5) makes sure that the infinite sum converges, and it is identical to the stability condition for the IIR filter given by (2).

Setting the derivative of (5) to zero results in

(6)(1α)N1(2α)2=1

Note that the optimal α must be in the interval (0,1] because larger values of α result in an alternating impulse response (3), which cannot approximate the constant impulse repsonse of the FIR moving average filter.

Taking the square root of (6) and introducing β=1α, we obtain

(7)β(N+1)/2+β(N1)/21=0

This equation cannot be solved analytically for β, but it can be solved for N:

(8)N=2log(1+β)log(β),β0

Equation (8) can be used to double-check a numerical solution of (7); it must return the specified value of N.

Equation (7) can be solved with a few lines of (Matlab/Octave) code:

N = 50;     % desired filter length of FIR moving average filter

if ( N == 1 )    % no iteration for trivial case
    b = 0;
else
    % Newton iteration
    b = 1;       % starting value
    Nit = 7;
    n = (N+1)/2;
    for k = 1:Nit,
        f = b^n + b^(n-1) -1;
        fp = n*b^(n-1) + (n-1)*b^(n-2);
        b = b - f/fp;
    end

    % check result
    N0 = -2*log(1+b)/log(b) + 1     % must equal N
end

a = 1 - b;

Below is a table with the optimal values of α for a range of filter lengths N:

   N     alpha

   1   1.0000e+00
   2   5.3443e-01
   3   3.8197e-01
   4   2.9839e-01
   5   2.4512e-01
   6   2.0809e-01
   7   1.8083e-01
   8   1.5990e-01
   9   1.4333e-01
  10   1.2987e-01
  20   6.7023e-02
  30   4.5175e-02
  40   3.4071e-02
  50   2.7349e-02
  60   2.2842e-02
  70   1.9611e-02
  80   1.7180e-02
  90   1.5286e-02
 100   1.3768e-02
 200   6.9076e-03 
 300   4.6103e-03
 400   3.4597e-03
 500   2.7688e-03
 600   2.3078e-03
 700   1.9785e-03
 800   1.7314e-03
 900   1.5391e-03
1000   1.3853e-03
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.