summaryrefslogtreecommitdiff
path: root/modules/signal_processing/macros/eqfir.sci
blob: d26a9132c49bd97ecb18bceea91d3f0b369c725f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - 1988 - C. Bunks
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt

function [hn]=eqfir(nf,bedge,des,wate)
    //<hn>=eqfir(nf,bedge,des,wate)
    //Minimax approximation of multi-band, linear phase, FIR filter
    //  nf    :Number of output filter points desired
    //  bedge :Mx2 matrix giving a pair of edges for each band
    //  des   :M-vector giving desired magnitude for each band
    //  wate  :M-vector giving relative weight of error in each band
    //  hn    :Output of linear-phase FIR filter coefficients
    //!

    //get number of cosines

    nc=int(nf/2);
    if nf-2*nc<>0 then,
        flag=0;
        nc=nc+1;
    else,
        flag=1;
    end,

    //make frequency grid, desired function, and weight function

    [nb,c2]=size(bedge);
    ngp=nc*16;
    b1=bedge(:,1);
    b2=bedge(:,2);
    sb=sum(b2-b1);
    delf=sb/ngp;
    bp=round((b2-b1)/delf);
    bsum=0;
    for k=1:nb,
        bpk=bp(k);
        et=b2(k)-b1(k);
        fg(bsum+1:bsum+bpk)=b1(k)*ones(1:bpk)+(0:bpk-1)*et/(bpk-1);
        ds(bsum+1:bsum+bpk)=des(k)*ones(1:bpk);
        wt(bsum+1:bsum+bpk)=wate(k)*ones(1:bpk);
        bsum=bsum+bpk;
    end,

    //adjust values of ds and wt if filter is of even length

    if flag==1 then,
        fgs=max(size(fg));
        if fg(fgs)>.5-%eps then,
            fg=fg(1:fgs-1);
            ds=ds(1:fgs-1);
            wt=wt(1:fgs-1);
        end,
        cf=cos(%pi*fg);
        ds=ds./cf;
        wt=wt.*cf;
    end,

    //call remez

    [an]=remezb(nc,fg,ds,wt);

    //obtain other half of filter coefficients (by symmetry)

    if flag==1 then,
        hn(1)=.25*an(nc);
        hn(2:nc-1)=.25*(an(nc:-1:3)+an(nc-1:-1:2));
        hn(nc)=.5*an(1)+.25*an(2);
        hn(nc+1:2*nc)=hn(nc:-1:1);
    else,
        hn=an(nc:-1:2)/2;
        hn(nc)=an(1);
        hn(nc+1:2*nc-1)=hn(nc-1:-1:1);
    end,
endfunction