Apply frequency band transformation from lowpass to bandpass in the z-domain.
Example 1
Elliptic bandpass filter design. The cutoff frequency of a lowpass analog filter prototype transformed in to z domain is obtained with BilinearUnwarp method. The analog prototype filter has a normalized cutoff frequency at 1 rad/sec.
uses MtxExpr, Math387, MtxVec, MtxVecTee, MtxVecEdit,
LinearSystems, IirFilters, SignalUtils;
procedure TForm1.Button1Click(Sender: TObject);
var z,p,num,den,Response: Vector;
Order: integer;
k,Wc,w1,w2,BW: TSample;
begin
Order := 4;//design a fourth order filter.
EllipticAnalog(Order,0.1,30,z,p,k); //design analog protype
Bilinear(z,p,k,2); //bilinear with sampling frequency 2Hz.
w1 := 0.2; //start of the passband at 0.2Hz.
w2 := 0.5; //stop of the passband at 0.5Hz.
Wc := Sqrt(w1*w2); //center frequency of the passband
BW := w2-w1; //passband width
// LowpassToBandpassZ(z,p,k,Wc,BW,BilinearUnwarp(1));
// ZeroPoleToTransferFun(num,den,z,p,k);
//Alternative:
// ...
ZeroPoleToTransferFun(num,den,z,p,k);
LowpassToBandpassZ(num,den,Wc,BW,BilinearUnwarp(1));
FrequencyResponse(num,den,Response,64);
DrawIt(Response);
end;
#include "MtxVecCPP.h" //MtxVecCPP.cpp must be included in the project
#include "MtxVecEdit.hpp"
#include "MtxVecTee.hpp"
#include "SignalUtils.hpp"
#include "IirFilters.hpp"
#include "LinearSystems.hpp"
void __fastcall TForm1::BitBtn1Click(TObject *Sender)
{
Vector z,p,num,den,Response;
int Order;
double k,Wc,w1,w2,BW;
Order = 4;//design a fourth order filter.
EllipticAnalog(Order,0.1,30,z,p,k); //design analog protype
Bilinear(z,p,k,2); //bilinear with sampling frequency 2Hz.
w1 = 0.2; //start of the passband at 0.2Hz.
w2 = 0.5; //stop of the passband at 0.5Hz.
Wc = Sqrt(w1*w2); //center frequency of the passband
BW = w2-w1; //passband width
// LowpassToBandpassZ(z,p,k,Wc,BW,BilinearUnwarp(1));
// ZeroPoleToTransferFun(num,den,z,p,k);
//Alternative:
// ...
ZeroPoleToTransferFun(num,den,z,p,k);
LowpassToBandPassZ(num,den,Wc,BW,BilinearUnwarp(1));
FrequencyResponse(num,den,Response,64);
DrawIt(Response);
}
using Dew.Math;
using Dew.Math.Editors;
using Dew.Math.Units;
using Dew.Signal;
using Dew.Signal.Units;
using Dew.Math.Tee;
using Dew.Signal.Tee;
private void button1_Click(object sender, EventArgs e)
{
Vector z = new Vector(0);
Vector p = new Vector(0);
Vector num = new Vector(0);
Vector den = new Vector(0);
Vector Response = new Vector(0);
double k, Wc;
double FS = 2;
int Order = 4; //design a fourth order filter.
IIRFilters.EllipticAnalog(Order,0.1,30, z, p, out k); //design analog protype
LinearSystems.Bilinear(z, p, ref k, FS,true);
double w1 = 0.2; //start of the passband at 0.2Hz.
double w2 = 0.5; //stop of the passband at 0.5Hz.
Wc = Math.Sqrt(w1*w2); //center frequency of the passband
double BW = w2-w1; //passband width
LinearSystems.LowpassToBandPassZ(z, p, ref k, Wc, BW, LinearSystems.BilinearUnwarp(1,FS));
LinearSystems.ZeroPoleToTransferFun(num,den, z, p, k);
SignalUtils.FrequencyResponse(num, den, Response, 64, false, TSignalWindowType.wtRectangular, 0); //zero padding set to 64
//Alternative:
// ...
// LinearSystems.ZeroPoleToTransferFun(num,den, z, p, k);
// LinearSystems.LowpassToBandPassZ(num,den, Wc, BW, LinearSystems.BilinearUnwarp(1,FS));
TeeChart.DrawIt(Response, "Frequency response", false);
//TeeChart.DrawIt(20 * MtxExpr.Log10(MtxExpr.Abs(Response)), "Magnitude", false);
//TeeChart.DrawIt(MtxExpr.PhaseSpectrum(Response) * (180 / Math.PI), "Phase", false);
}
Declaration
Procedure LowpassToBandPassZ(num, den: TVec; Freq, BW: TSample; PrototypeFreq: TSample);
Description
Freq is the center frequency of the passband with width BW of the new filter. The function returns modified num and den. PrototypeFreq is the cutoff frequency of the prototype lowpass filter after it has been mapped to z-domain. Freq, BW and PrototypeFreq must be between 0 and 1 (Sampling frequency = 2). The transformation is defined with the following mapping ([1] p. 260 and [2] p. 434, [3] p. 352):
-a2 + a1*z^(-1) - z^(-2)
z^(-1) ---> ---------------------------
1 - a1*z^(-1) + a2*z^(-2)
cos((w2 + w1)/2)
Beta = -------------------
cos((w2 - w1)/2)
k = cot((w2-w1)/2)*tan(wc/2)
a1 = 2*beta*k1/(k1+1)
a2 = (k1-1)/(k1+1)
wc - old cutoff frequency
w2 - desired upper cutoff frequency
w1 - desired lower cutoff frequency