Remez Octave Signal vs Scipi.signal

[Thomas D. Dean]

Here is a test I tried.  The lpf values are from scipy.signal.remez.  It 
takes several minutes to run.  (150% cpu 4.2 GHz i7)

Tom Dean

## Compare the remez.cc output to scipy.signal.remez
1; ## not a function file
##
## from sdr_fm_decode.py
##
## Design filter taps=64 200000 292500 570000
## lpf = signal.remez(n_taps, [0, f_bw, 
f_bw+(Fs/2-f_bw)/4,Fs/2],[1,0],Hz=Fs)
f_bw = 200000; ## bandwidth
Fs = 1140000;  ## sample rate
n_taps = 64;

lpf =[ -5.94491898e-05, -6.29398695e-05,  1.36535572e-04,  2.23759914e-04,
       -1.49209793e-04, -5.66615829e-04, -3.31926564e-05   1.02819068e-03,
        6.43169323e-04, -1.41717579e-03, -1.84240903e-03,  1.30670204e-03,
        3.61493083e-03, -1.25323608e-04, -5.56377435e-03, -2.69078248e-03,
        6.82099667e-03,  7.39922235e-03, -6.06610456e-03, -1.36412904e-02,
        1.72606880e-03,  2.01803166e-02,  7.74640578e-03, -2.47593374e-02,
       -2.37303856e-02,  2.38912118e-02,  4.80521101e-02, -1.15316628e-02,
       -8.70066265e-02, -3.09559112e-02,  1.88176232e-01,  3.99235995e-01,
        3.99235995e-01,  1.88176232e-01, -3.09559112e-02, -8.70066265e-02,
       -1.15316628e-02,  4.80521101e-02,  2.38912118e-02, -2.37303856e-02,
       -2.47593374e-02,  7.74640578e-03,  2.01803166e-02,  1.72606880e-03,
       -1.36412904e-02, -6.06610456e-03,  7.39922235e-03,  6.82099667e-03,
       -2.69078248e-03, -5.56377435e-03, -1.25323608e-04,  3.61493083e-03,
        1.30670204e-03, -1.84240903e-03, -1.41717579e-03,  6.43169323e-04,
        1.02819068e-03, -3.31926564e-05  -5.66615829e-04, -1.49209793e-04,
        2.23759914e-04,  1.36535572e-04, -6.29398695e-05  -5.94491898e-05];
## change to a column vector
lpf = reshape(lpf,numel(lpf),1);

## The band edges need to be changed to the range 0..1 and there needs
## to be one amplitude for each band edge.

clear remez
filt = remez(n_taps,
             [0, f_bw, f_bw+(Fs/2-f_bw)/4,Fs/2]/(Fs/2),
             [1,1,0,0]);

dat = [ones(n_taps/2,1);zeros(n_taps/2,1)];
ifilt = ifft(dat);

figure()
hold on
subplot(3,1,1)
plot(lpf)
title('Python scipy.signal.remez')
subplot(3,1,2)
plot(filt)
title('Octave remez.cc')
subplot(3,1,3)
plot(abs(ifilt))
title('ifft')
hold off

num_samp = 8192000;
x2 = rand(num_samp,1);

## apply filter
x3_lpf=filtfilt(lpf,1,x2);
x3_filt=filtfilt(filt,1,x2);

## decimate
dec_rate = floor(Fs / f_bw);
x4_lpf = x3_lpf(1:dec_rate:end);
x4_filt = x3_filt(1:dec_rate:end);

## new sample rate
Fs_y = Fs / dec_rate;

figure()
hold on
subplot(3,1,1)
pwelch(x2,4,0,2048,Fs)
title('pwelch x2')
subplot(3,1,2)
pwelch(x4_lpf,4,0,2048,Fs_y)
title('pwelch x4_lpf')
subplot(3,1,3)
pwelch(x4_filt,4,0,2048,Fs_y)
title('pwelch x4_filt')
hold off