Remez Octave Signal vs Scipi.signal

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

Remez Octave Signal vs Scipi.signal

Thomas D. Dean-2
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

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Remez Octave Signal vs Scipi.signal

Dmitri A. Sergatskov
On Sun, Sep 3, 2017 at 10:49 PM, Thomas D. Dean <[hidden email]> wrote:

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



​I am not a digital filter guy, but
pwelch(x4_lpf,4,0,2048,Fs_y)

​does not make any sense to me.


​Try
​pwelch(x4_lpf,[],[],2048,Fs_y, "loglog", "linear")

(you pick up 4 data-point window, pad it with 0 to 2048, do fft on that,
make a zero overlap -- what is the point of pwelch then?)

Also there is a decimate() function in signal pkg.

​Dmitri.
--


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Remez Octave Signal vs Scipi.signal

Thomas D. Dean-2
On 09/03/2017 09:28 PM, Dmitri A. Sergatskov wrote:

> On Sun, Sep 3, 2017 at 10:49 PM, Thomas D. Dean <[hidden email]
> <mailto:[hidden email]>> wrote:
>
>
>     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
>
>
>
> ​I am not a digital filter guy, but
> ​
> pwelch(x4_lpf,4,0,2048,Fs_y)
>
> ​does not make any sense to me.
>
>
> ​Try
> ​pwelch(x4_lpf,[],[],2048,Fs_y, "loglog", "linear")
>
> (you pick up 4 data-point window, pad it with 0 to 2048, do fft on that,
> make a zero overlap -- what is the point of pwelch then?)
>
> Also there is a decimate() function in signal pkg.

Thanks,  I know less about DSP than almost everyone.  I did the math for
DSP 40 years ago at Cal.

The difficulty I had was in the values returned by the
Octave:signal:remez function as compared to the scipy.signal.remez function.

I don't remember why one cannot do an fft, zero all the data that does
not fit the filter criteria, and do an ifft.  The math must back this
up, but, I don't remember it...  I am looking for a book that says, "if
you want that result, do this."

Tom Dean

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave