diff --git a/applications/lofar2/model/rtdsp/firfilter.py b/applications/lofar2/model/rtdsp/firfilter.py index cbfb0a39150e11b474427ad7cb9fa478d5c1c2d2..beea36c2d2f5758d9f5973ef6d46a528a4a88092 100644 --- a/applications/lofar2/model/rtdsp/firfilter.py +++ b/applications/lofar2/model/rtdsp/firfilter.py @@ -156,65 +156,10 @@ def ideal_low_pass_filter(Npoints, Npass, bandEdgeGain=1.0): return h, f, HF -def prototype_fir_low_pass_filter(method='firls', - Npoints=1024, Ntaps=16, Ncoefs=1024*16, - hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, kaiserBeta=1, fs=1.0): - """Derive FIR coefficients for prototype low pass filter - - Use method 'firls' or 'remez'. - - Default use LPF specification for LOFAR subband filter. For subband filter - BWbin = fs / Npoints and the number of subbands is Nsub = Npoints / 2. - - For multirate upsampling or or downsampling filter the Npoints is the - multirate factor Nup in upsampling or Ndown in downsampling. The fpass = - BWbin / 2 = fNyquist / Npoints, where fNyquist = fs / 2. - - Use hpFactor < 1.0 to increase attentuation at fpass, to reduce disturbance - due to aliasing. - - The number of prototype FIR filter coefficients Ncoefs <= Npoints * Ntaps. - For symmetrical-coefficient FIR filters the group delay is - (Ncoefs - 1) / 2, so choose Ncoefs is odd to have integers number of - samples output delay. - - Normalize DC gain to 1.0. - - Input: - - method: FIR filter design method 'firls' or 'remez' - - Npoints: Number of points of the DFT in the filterbank, or multirate - factor Nup in upsampling or Ndown in downsampling - - Ntaps: Number of taps per polyphase. - - Ncoefs: Actual number of prototype FIR filter coefficients, any extra - coefs are zero - - hpFactor : Half power bandwidth of the filter relative to BWbin - - transitionFactor: transition bandwidth factor relative to fpass - - stopRippleFactor: stopband ripple factor relative to pass band ripple - - kaiserBeta: When kaiserBeta > 0 then additionally apply a Kaiser window on FIR - coefficients - - fs: sample frequency, for logging - Return: - - h: FIR coefficients for requested Ncoefs - """ - # LPF specification for subband filter - BWbin = fs / Npoints # bandwidth of one bin (is subband frequency) - # Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0 - BWpass = hpFactor * BWbin - fpass = BWpass / 2 # pass band frequency of bin at DC: -fpass to +fpass - transitionBW = transitionFactor * fpass # transition bandwidth - fcutoff = 0 # not used, use hpFactor instead - cutoffGain = 0.5 - fstop = fpass + transitionBW # stop band frequency - rippleWeights = [1, stopRippleFactor] - - # Design subband filter - h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, kaiserBeta, fs) - return h - - def design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, fcutoff=0, cutoffGain=0.5, rippleWeights=[1, 1], - kaiserBeta=0, fs=1.0): + kaiserBeta=0, fs=1.0, + cutoffBW=0): """Derive FIR coefficients for prototype low pass filter Use method 'firls' or 'remez', fs = 1.0 @@ -225,6 +170,17 @@ def design_fir_low_pass_filter(method, Normalize DC gain to 1.0. + Remark: + . The actual LPF gain at fcutoff will differ from specified cutoffGain by about 1%. + + Workaround: + . Method 'firls' does not support zero BW at fcutoff. With single fcutoff point error occurs: 'bands must be + monotonically nondecreasing and have width > 0'. Work around: use a very narrow band cutoffBW at fcutoff. + The cutoffBW band has to be much smaller than the transition band transistionBW, to avoid a flattened gain + around fcutoff, but also not too small, to because then the cutoffGain at fcutoff is not achieved. + Default force cutoffBW = transistionBW / 10 for 'firls' when input cutoffBW=0. For 'remez' method the + cutoffBW=0 effectively, without work around, because 'remez' does support single fcutoff value specification. + Input: - method: FIR filter design method 'firls' or 'remez' - Ncoefs: Actual number of prototype FIR filter coefficients, any extra @@ -237,11 +193,13 @@ def design_fir_low_pass_filter(method, - rippleWeights: relative ripple factors for pass band, optional fcutoff, and stop band - kaiserBeta: When kaiserBeta > 0, then additionally apply a Kaiser window on FIR coefficients + - cutoffBW: workaround bandwidth at fcutoff for firls, because firls does not support zero BW at fcutoff. Return: - h: FIR coefficients for requested Ncoefs """ fNyquist = fs / 2 - fsub = fpass * 2 # double sided LPF band width, when used as prototype for a BPF + lpBW = fpass + fstop # low pass bandwidth is twice center of transition band + transistionBW = fstop - fpass # Initial FIR filter length N # . Use interpolation of factor Q shorter filter to ensure the FIR filter design converges @@ -268,6 +226,7 @@ def design_fir_low_pass_filter(method, f_pb = fpass * Q # pass band frequency f_co = fcutoff * Q # optional cutoff frequency in transition band f_sb = fstop * Q # stop band frequency + f_tb = transistionBW * Q # Calculate initial FIR filter for length N if method == 'firls': @@ -275,7 +234,11 @@ def design_fir_low_pass_filter(method, hFir = signal.firls(N, [0, f_pb, f_sb, fNyquist], [1, 1, 0, 0], rippleWeights, fs=fs) else: - hFir = signal.firls(N, [0, f_pb, f_co, f_co, f_sb, fNyquist], + if cutoffBW == 0: + b_co = f_tb / 10 # default + else: + b_co = cutoffBW * Q # from input argument + hFir = signal.firls(N, [0, f_pb, f_co - b_co / 2, f_co + b_co / 2, f_sb, fNyquist], [1, 1, cutoffGain, cutoffGain, 0, 0], rippleWeights, fs=fs) if method == 'remez': if not fcutoff: @@ -305,18 +268,74 @@ def design_fir_low_pass_filter(method, h = fourier_interpolate(HFfir, Ncoefs) print('> design_fir_low_pass_filter()') - print('. Method = %s' % method) - print('. Q = %f' % Q) - print('. fsub = fpass * 2 = %f' % fsub) - print('. fpass = %f' % fpass) + print('. Method = %s' % method) + print('. Q = %f' % Q) + print('. fpass = %f' % fpass) + print('. fstop = %f' % fstop) + print('. lpBW = %f' % lpBW) if fcutoff: - print('. fcutoff = %f' % fcutoff) - print('. cutoffGain = %f' % cutoffGain) - print('. fstop = %f' % fstop) - print('. fNyquist = %f' % fNyquist) - print('. fs = %f' % fs) - print('. Ncoefs = %d' % len(h)) - print('. DC sum = %f' % np.sum(h)) - print('. Symmetrical coefs = %s' % is_symmetrical(h)) + print('. fcutoff = %f' % fcutoff) + print('. cutoffGain = %f' % cutoffGain) + print('. fNyquist = %f' % fNyquist) + print('. fs = %f' % fs) + print('. Ncoefs = %d' % len(h)) + print('. DC sum = %f' % np.sum(h)) + print('. Symmetrical coefs = %s' % is_symmetrical(h)) print('') return h + + +def prototype_fir_low_pass_filter(method='firls', + Npoints=1024, Ntaps=16, Ncoefs=1024*16, + hpFactor=0.9, transitionFactor=0.4, stopRippleFactor=1000000, kaiserBeta=0, fs=1.0): + """Derive FIR coefficients for prototype low pass filter + + Use method 'firls' or 'remez'. + + Default use LPF specification for LOFAR subband filter. For subband filter + BWbin = fs / Npoints and the number of subbands is Nsub = Npoints / 2. + + For multirate upsampling or or downsampling filter the Npoints is the + multirate factor Nup in upsampling or Ndown in downsampling. The fpass = + BWbin / 2 = fNyquist / Npoints, where fNyquist = fs / 2. + + Use hpFactor < 1.0 to increase attentuation at fpass, to reduce disturbance + due to aliasing. + + The number of prototype FIR filter coefficients Ncoefs <= Npoints * Ntaps. + For symmetrical-coefficient FIR filters the group delay is + (Ncoefs - 1) / 2, so choose Ncoefs is odd to have integers number of + samples output delay. + + Normalize DC gain to 1.0. + + Input: + - method: FIR filter design method 'firls' or 'remez' + - Npoints: Number of points of the DFT in the filterbank, or multirate + factor Nup in upsampling or Ndown in downsampling + - Ntaps: Number of taps per polyphase. + - Ncoefs: Actual number of prototype FIR filter coefficients, any extra + coefs are zero + - hpFactor : Half power bandwidth of the filter relative to BWbin + - transitionFactor: transition bandwidth factor relative to fpass + - stopRippleFactor: stopband ripple factor relative to pass band ripple + - kaiserBeta: When kaiserBeta > 0 then additionally apply a Kaiser window on FIR + coefficients + - fs: sample frequency, for logging + Return: + - h: FIR coefficients for requested Ncoefs + """ + # LPF specification for subband filter + BWbin = fs / Npoints # bandwidth of one bin (is subband frequency) + # Use half power bandwidth factor to tune half power cutoff frequency of LPF, default 1.0 + BWpass = hpFactor * BWbin + fpass = BWpass / 2 # pass band frequency of bin at DC: -fpass to +fpass + transitionBW = transitionFactor * fpass # transition bandwidth + fcutoff = 0 # not used, use hpFactor instead + cutoffGain = 0.5 + fstop = fpass + transitionBW # stop band frequency + rippleWeights = [1, stopRippleFactor] + + # Design subband filter + h = design_fir_low_pass_filter(method, Ncoefs, fpass, fstop, fcutoff, cutoffGain, rippleWeights, kaiserBeta, fs) + return h