FIR, band pass, band stop, high pass

% this is BANDPASS

close all; clear all; clc;

fm = 16e3;
fp1 = 0.7e3;
fp2 = 3.5e3;
fs1 = 0.4e3;
fs2 = 4.5e3;

Ap = 1;
As = 45;

Fp1 = fp1/fm;
Fp2 = fp2/fm;
Fs1 = fs1/fm;
Fs2 = fs2/fm;

Flpp = (Fp2 - Fp1)/2;
Flps = (Fs2 - Fs1)/2;

F0 = (Fp2+Fp1)/2;

N = 9;
Fc = Flpp;
dFc = (Flps-Flpp)/50;
run = 1;
while run ==1
    n = 0:N-1;
    w = hanning(N);
    hlp = 2*Fc*sinc(2*Fc*(n-(N-1)/2)).*w';
    At = freqz(hlp,1,[Flpp Flps],1);
    Atdb = 20*log10(abs(At));

    Apass = Atdb(1);   % passband edge
    Astop = Atdb(2);   % stopband edge

    if (Apass < -Ap) || (Astop > -As)
        if Fc+dFc >Flps
            N = N+2;
            Fc = Flpp;
        else
            Fc = Fc + dFc;
        end
    else
        run = 0;
    end
end

h = 2*cos(2*pi*F0*(n-(N-1)/2)).*hlp;

figure();
freqz(h,1,5000,fm);
figure();
grpdelay(h,1,5000,fm);
figure();
impz(h,1,40);


% This is BANDSTOP

close all; clear all; clc;

fm = 16e3;
fp1 = 0.7e3;
fp2 = 3.5e3;
fs1 = 0.4e3;
fs2 = 4.5e3;

Ap = 1;
As = 45;

Fp1 = fp1/fm;
Fp2 = fp2/fm;
Fs1 = fs1/fm;
Fs2 = fs2/fm;

Flpp = (Fp2 - Fp1)/2;
Flps = (Fs2 - Fs1)/2;

F0 = (Fp2+Fp1)/2;

N = 9;
Fc = Flpp;
dFc = (Flps-Flpp)/50;
run = 1;

while run == 1
    n = 0:N-1;
    w = hanning(N);
    hlp = 2*Fc*sinc(2*Fc*(n-(N-1)/2)).*w';

    At = freqz(hlp,1,[Flpp Flps],1);
    Atdb = 20*log10(abs(At));

    Apass = Atdb(1);
    Astop = Atdb(2);

    if (Apass < -Ap) || (Astop > -As)
        if Fc+dFc > Flps
            N = N+2;
            Fc = Flpp;
        else
            Fc = Fc + dFc;
        end
    else
        run = 0;
    end
end

hbp = 2*cos(2*pi*F0*(n-(N-1)/2)).*hlp;
h = -hbp;
h((N+1)/2) = h((N+1)/2) + 1;

figure; freqz(h,1,5000,fm);
figure; grpdelay(h,1,5000,fm);
figure; impz(h,1,40);


% this is high pass

close all; clear all; clc;

fm = 16e3;
fp = 3.5e3;
fs = 0.7e3;

Ap = 1;
As = 45;

Fp = fp/fm;
Fs = fs/fm;

Flpp = Fs;
Flps = Fp;

N = 9;
Fc = Flpp;
dFc = (Flps-Flpp)/50;
run = 1;

while run == 1
    n = 0:N-1;
    w = hanning(N);
    hlp = 2*Fc*sinc(2*Fc*(n-(N-1)/2)).*w';

    At = freqz(hlp,1,[Flpp Flps],1);
    Atdb = 20*log10(abs(At));

    Apass = Atdb(1);
    Astop = Atdb(2);

    if (Apass < -Ap) || (Astop > -As)
        if Fc+dFc > Flps
            N = N+2;
            Fc = Flpp;
        else
            Fc = Fc + dFc;
        end
    else
        run = 0;
    end
end

h = -hlp;
h((N+1)/2) = h((N+1)/2) + 1;

figure; freqz(h,1,5000,fm);
figure; grpdelay(h,1,5000,fm);
figure; impz(h,1,40);

Answer 1

[Mathieu NOE]

hello

your code is correct but gives some warnings :

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.189727e-64.

> In bilinear (line 169)

In Cheb1 (line 22)

NB that there is a much more direct way to design digital Cheb filters - I let you overlay the results (vs your code) :

rp=1;   % max passband ripple (in dB)
rs=40;  % min stopband attenuation (in dB)
fp=20e3;   % passband freq
fs=16e3;   % stopband freq
fm=64e3;    % sampling freq

% normalize freqs vs Nyquist
Nwp=2*fp/fm;
Nws=2*fs/fm;

[N,wc]=cheb1ord(Nwp,Nws,rp,rs);
[Bz,Az]=cheby1(N,rp,wc,'high');     % high pass

figure();
freqz(Bz,Az,1024,fm);
title(sprintf('N=%d cheby-type-1 Hp IIR filter',N));

●

figure();
grpdelay(Bz, Az, 1024, Az);

figure();
impz(Bz,Az,50);