头文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* advmath.h
*
* Some advance math functions used in IIR filter design.
*
* Adapted form "Analog and Digital Filter Design", Les Thede, 2004.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef ADVMATH_H
#define ADVMATH_H
#include <cmath>
#include <complex>
namespace splab
{
inline double acosh( double );
inline double asinh( double );
double arcsc( double, double );
void ellipticFun( double, double, double*, double*, double* );
double ellipticIntegral( double );
void quadradicRoot( complex<double>, complex<double>, complex<double>,
complex<double>&, complex<double>& );
#include <advmath-impl.h>
}
// namespace splab
#endif
// ADVMATH_H/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* iir.h
*
* Infinite Impulse Response.
*
* This class is designed for designing the IIR filter using bilinear
* transform method. The unit of frequency and gain parameters are "Hz" and
* "dB", respectively.
*
* The valid filter types are:
* lowpass highpass bandpass bandstop
* and the valid approximation method are:
* Butterworth Chebyshev Inverse Chebyshev Elliptic
*
* The filter's coefficients consist of quadratic terms( the linear term can
* be viewed as: a + b*z^-1 + 0*z^-2 ). The numerator and denominator
* coefficients are stored in "aCoefs" and "bCoefs", respectively. The kth
* quadratic term is as follow:
*
* aCoefs[3k] + aCoefs[3k+1]*z^-1 + aCoefs[3k+2]*z^-2
* -------------------------------------------------------
* bCoefs[3k] + bCoefs[3k+1]*z^-1 + bCoefs[3k+2]*z^-2
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef IIR_H
#define IIR_H
#include <iomanip>
#include <cstdio>
#include <dfd.h>
#include <advmath.h>
namespace splab
{
using std::ios;
using std::setw;
using std::setiosflags;
using std::setprecision;
class IIR : public DFD
{
public:
IIR( const string &select, const string &appMethod );
~IIR();
void design();
void dispInfo() const;
Vector<double> getNumCoefs() const;
Vector<double> getDenCoefs() const;
private:
void freqWrap();
bool orderEst();
void normCoefs();
void unnormCoefs();
void bilinTran();
void freqUnwrap();
void butterCoefs();
void chebyCoefs();
void invChebyCoefs();
void elliptCoefs();
void unnormLP( double freq );
void unnormHP( double freq );
void unnormBP( double BW, double Wo );
void unnormBS( double BW, double Wo );
double frqeResp( double freq );
void calcGain();
string approx; // approximation method
double gain; // overall gain
Vector<double> aCoefs, // numerator
bCoefs, // denominator
edgeGain; // gain at edge frquency
};
// class IIR
#include <iir-impl.h>
}
// namespace splab
#endif
// IIR_H实现文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* advmath-impl.h
*
* Implementationfor advance math functions.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
/**
* Calcutes the inverse hyperbolic cosine of x.
*/
inline double acosh( double x )
{
return log( x + sqrt(x*x-1) );
}
/**
* Calculates the inverse hyperbolic sine of x.
*/
inline double asinh( double x )
{
return log( x + sqrt(x*x + 1) );
}
/**
* Calculates Jacobian elliptic arctan function using
* arithmetic-geometric mean method.
*/
double arcsc( double u, double k )
{
int i, L;
double A, B, BT, Y;
A = 1;
B = k;
Y = 1.0 / u;
L = 0;
for( i=0; i<MAXTERM; ++i )
{
BT = A * B;
A = A + B;
B = 2 * sqrt(BT);
Y = Y - BT / Y;
if( Y == 0 )
Y = sqrt(BT) * EPS;
if( abs(A-B) < (A*EPS) )
break;
L = 2 * L;
if( Y < 0 )
L++;
}
if( Y < 0 )
L++;
return ( (atan(A/Y) + PI * L) / A );
}
/**
* Calculates Jacobian elliptic functions using
* arithmetic-geometric mean method.
*/
void ellipticFun( double u, double k,
double *sn, double *cn, double *dn )
{
int i, imax;
double A[MAXTERM], B[MAXTERM],
C[MAXTERM], P[MAXTERM];
k = k * k;
A[0] = 1;
B[0] = sqrt(1-k);
C[0] = sqrt(k);
for( i=1; i<MAXTERM; ++i )
{
A[i] = (A[i-1] + B[i-1]) / 2;
B[i] = sqrt( A[i-1] * B[i-1] );
C[i] = (A[i-1] - B[i-1]) / 2;
if( C[i] < EPS )
break;
}
if( i == MAXTERM )
imax = i - 1;
else
imax = i;
P[imax] = pow(2.0,imax) * A[imax] * u;
for( i=imax; i>0; i-- )
P[i-1] = (asin(C[i]*sin(P[i])/A[i]) + P[i]) / 2;
*sn = sin(P[0]);
*cn = cos(P[0]);
*dn = sqrt( 1 - k*(*sn)*(*sn) );
}
/**
* Calculates complete elliptic integral using
* arithmetic-geometric mean method.
*/
double ellipticIntegral( double k )
{
int i;
double A[MAXTERM], B[MAXTERM], C[MAXTERM];
k = k*k;
A[0] = 1;
B[0] = sqrt(1-k);
C[0] = sqrt(k);
for( i=1; i<MAXTERM; ++i )
{
A[i] = (A[i-1] + B[i-1]) / 2;
B[i] = sqrt(A[i-1]*B[i-1]);
C[i] = (A[i-1] - B[i-1]) / 2;
if( C[i] < EPS )
break;
}
return PI / (2*A[i]);
}
/**
* Factoring quadratic equation with complex coefficients.
* Equation form is a*x^2 + b*x + c, solutions will be placed
* in complex numbers pointers d and e.
*/
void quadradicRoot( complex<double> a, complex<double> b, complex<double> c,
complex<double> &x1, complex<double> &x2 )
{
complex<double> a2,
ac4,
sq;
a2 = 2.0 * a;
sq = sqrt( b*b - 4.0*a*c );
x1 = ( -b + sq ) / a2;
x2 = ( -b - sq ) / a2;
}/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* iir-impl.h
*
* Implementation for IIR class.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
/**
* constructors and destructor
*/
IIR::IIR( const string &select, const string &app )
: DFD( select ), approx(app)
{
bool cond = ( filtType=="lowpass" || filtType=="highpass" ||
filtType=="bandpass" || filtType=="bandstop" );
assert(cond);
cond = ( approx=="Butterworth" || approx=="Chebyshev" ||
approx=="InvChebyshev" || approx=="Elliptic" );
assert(cond);
}
IIR::~IIR()
{
}
/**
* Calculates the digital FIR coefficients.
*/
void IIR::design()
{
freqWrap();
if( orderEst() )
{
normCoefs();
unnormCoefs();
bilinTran();
freqUnwrap();
calcGain();
}
else
cerr << "The filter's order exceeds exceed the range!" << endl;
}
/**
* Displays the design result of filter.
*/
void IIR::dispInfo() const
{
int i, j;
cout << endl;
cout << "\t\t Filter selectivity : " << filtType << endl;
cout << "\t\t Approximation method : " << approx << endl;
cout << "\t\t Filter order : " << order << endl;
cout << "\t\t Overall gain : " << gain << endl;
cout << "\t\t Sampling Frequency (Hz) : " << fsamp << endl;
// gains and edge frequency
if( filtType=="lowpass" )
{
cout << "\t\t Passband frequency (Hz) : "
<< wpass1/TWOPI << endl;
cout << "\t\t Passband gain (dB) : "
<< apass1 << endl;
cout << "\t\t Stopband frequency (Hz) : "
<< wstop1/TWOPI << endl;
cout << "\t\t Stopband gain (dB) : "
<< astop1 << endl;
}
else if( filtType=="highpass" )
{
cout << "\t\t Stopband frequency (Hz) : "
<< wstop1/TWOPI << endl;
cout << "\t\t Stopband gain (dB) : "
<< astop1 << endl;
cout << "\t\t Passband frequency (Hz) : "
<< wpass1/TWOPI << endl;
cout << "\t\t Passband gain (dB) : "
<< apass1 << endl;
}
else if( filtType=="bandpass" )
{
cout << "\t\t Lower stopband frequency (Hz) : "
<< wstop1/TWOPI << endl;
cout << "\t\t Lower stopband gain (dB) : "
<< astop1 << endl;
cout << "\t\t Lower passband frequency (Hz) : "
<< wpass1/TWOPI << endl;
cout << "\t\t Upper passband frequency (Hz) : "
<< wpass2/TWOPI << endl;
cout << "\t\t Passband gain (dB) : "
<< apass1 << endl;
cout << "\t\t Upper stopband frequency (Hz) : "
<< wstop2/TWOPI << endl;
cout << "\t\t Upper stopband gain (dB) : "
<< astop2 << endl;
}
else
{
cout << "\t\t Lower passband frequency (Hz) : "
<< wpass1/TWOPI << endl;
cout << "\t\t Lower passband gain (dB) : "
<< apass1 << endl;
cout << "\t\t Lower stopband frequency (Hz) : "
<< wstop1/TWOPI << endl;
cout << "\t\t Upper stopband frequency (Hz) : "
<< wstop2/TWOPI << endl;
cout << "\t\t Stopband gain (dB) : "
<< astop1 << endl;
cout << "\t\t Upper passband frequency (Hz) : "
<< wpass2/TWOPI << endl;
cout << "\t\t Upper passband gain (dB) : "
<< apass2 << endl;
}
// display coefficients
cout << endl << endl << " Numerator Coefficients ";
cout << " Denominator Coefficients" << endl << endl ;
cout << " N [ 1 z^-1 z^-2 ]";
cout << " [ 1 z^-1 z^-2 ]" << endl;
cout << " === =================================";
cout << " =================================" << endl ;
for( i=0; i<(order+1)/2; ++i )
{
j = i * 3;
cout << setw(3) << i+1 << " ";
cout << setiosflags(ios::fixed );
cout << setw(6) << setprecision(1) << aCoefs[j+0];
cout << setw(14) << setprecision(8) << aCoefs[j+1];
cout << setw(14) << setprecision(8) << aCoefs[j+2];
cout << setw(10) << setprecision(1) << bCoefs[j+0];
cout << setw(14) << setprecision(8) << bCoefs[j+1];
cout << setw(14) << setprecision(8) << bCoefs[j+2] << endl ;
}
// display edge frequency response
cout << endl << endl
<< "\t ==================== Edge Frequency Response";
cout << " ====================" << endl;
cout << setiosflags(ios::fixed);
if( edgeGain.size() == 2 )
{
cout << "\t Mag(fp) = " << edgeGain[0] << "(dB)";
cout << " Mag(fs) = " << edgeGain[1] << "(dB)" << endl;
}
else
{
cout << "\t Mag(fp1) = " << edgeGain[0] << "(dB)";
cout << " Mag(fp2) = " << edgeGain[1] << "(dB)" << endl;
cout << "\t Mag(fs1) = " << edgeGain[2] << "(dB)";
cout << " Mag(fs2) = " << edgeGain[3] << "(dB)" << endl;
}
}
/**
* Get the FIR filter's coefficients.
*/
inline Vector<double> IIR::getNumCoefs() const
{
return aCoefs;
}
inline Vector<double> IIR::getDenCoefs() const
{
return bCoefs;
}
/**
* Pre-warps analog frequency for IIR filtrs.
*/
void IIR::freqWrap( )
{
// Pre-warps critical frequency based on selectivity.
if( filtType=="lowpass" || filtType=="highpass" )
{
wpass1 = 2*fsamp * tan( wpass1/(2*fsamp) );
wstop1 = 2*fsamp * tan( wstop1/(2*fsamp) );
}
else
{
wpass1 = 2*fsamp * tan( wpass1/(2*fsamp) );
wpass2 = 2*fsamp * tan( wpass2/(2*fsamp) );
wstop1 = 2*fsamp * tan( wstop1/(2*fsamp) );
wstop2 = 2*fsamp * tan( wstop2/(2*fsamp) );
}
}
/**
* Calculate the order of digital IIR filter.
*/
bool IIR::orderEst()
{
double kernel, ratio,
n,
rt, kn,
CEIrt, CEIrtq,
CEIkn, CEIknq,
wp1, wp2, ws1, ws2;
// Simplify the variables.
wp1 = wpass1;
wp2 = wpass2;
ws1 = wstop1;
ws2 = wstop2;
// Calculate frequency ratio based on filter selectivity.
if( filtType == "lowpass" )
ratio = ws1 / wp1;
else if( filtType == "highpass" )
ratio = wp1 / ws1;
else if( filtType == "bandpass" )
{
if( ws1 > (wp1*wp2)/ws2 )
{
ws2 = (wp1*wp2) / ws1;
wstop2 = ws2;
}
else
{
ws1 = (wp1*wp2) / ws2;
wstop1 = ws1;
}
ratio = (ws2-ws1) / (wp2-wp1);
}
else
{
if( wp1 > (ws1*ws2)/wp2 )
{
wp2 = (ws1*ws2) / wp1;
wpass2 = wp2;
}
else
{
wp1 = (ws1*ws2) / wp2;
wpass1 = wp1;
}
ratio = (wp2-wp1) / (ws2-ws1);
}
// Calculate kernel which is used in all calculations.
kernel = ( pow(10.0,-0.1*astop1) - 1 )
/ ( pow(10.0,-0.1*apass1) - 1 );
// Calculate order and round up to integer value.
if( approx == "Butterworth" )
n = log10(kernel) / (2*log10(ratio));
else if( approx == "Chebyshev" || approx == "InvChebyshev" )
n = acosh(sqrt(kernel)) / acosh(ratio);
else
{
rt = 1 / ratio;
if( rt > 0.9999 )
return false;
kn = 1 / sqrt(kernel);
if( kn < 2e-8 )
return false;
CEIrt = ellipticIntegral(rt);
CEIrtq = ellipticIntegral(sqrt(1-rt*rt));
CEIkn = ellipticIntegral(kn);
CEIknq = ellipticIntegral(sqrt(1-kn*kn));
n = ( CEIrt*CEIknq ) / ( CEIrtq * CEIkn );
}
if( n > 200 )
return false;
order = int( std::ceil(n) );
return true;
}
/**
* Calculate normal filter coefficients.
*/
void IIR::normCoefs()
{
// There are 3 coefficients for each quadratic.
// First order factors are considered as quadratics.
int coefsNum = 3 * ( (order+1) / 2 );
aCoefs.resize(coefsNum);
bCoefs.resize(coefsNum);
// Calculate coefs based on approximation.
if( approx == "Butterworth" )
butterCoefs();
else if( approx == "Chebyshev" )
chebyCoefs();
else if( approx == "InvChebyshev" )
invChebyCoefs();
else
elliptCoefs();
}
/**
* Converts normal lowpass coefficients to unnormalized LP/HP/BP/BS.
*/
void IIR::unnormCoefs()
{
double freq, BW, Wo;
// Calculate frequency, Wo and BW based on approx method.
if( approx == "InvChebyshev" )
{
freq = wstop1;
Wo = sqrt( wstop1 * wstop2 );
BW = wstop2 - wstop1;
}
else
{
freq = wpass1;
Wo = sqrt( wpass1 * wpass2 );
BW = wpass2 - wpass1;
}
// Call unnormal function based on selectivity.
if( filtType == "lowpass" )
unnormLP( freq );
else if( filtType == "highpass" )
unnormHP( freq );
else if( filtType == "bandpass" )
unnormBP( BW, Wo );
else
unnormBS( BW, Wo );
}
/**
* Use bilinear transform to convert transfer function
* from s-domain to z-domain.
*/
void IIR::bilinTran()
{
int i, j, start, numQuads;
double f2, f4,
N0, N1, N2,
D0, D1, D2;
f2 = 2 * fsamp;
f4 = f2 * f2;
numQuads = ( order + 1 ) / 2;
// Handle first order factor if present start
// indicates if one quad factor handled.
start = 0;
if( order % 2 )
{
N0 = aCoefs[2] + aCoefs[1] * f2;
N1 = aCoefs[2] - aCoefs[1] * f2;
D0 = bCoefs[2] + bCoefs[1] * f2;
D1 = bCoefs[2] - bCoefs[1] * f2;
aCoefs[0] = 1.0;
aCoefs[1] = N1 / N0;
aCoefs[2] = 0.0;
bCoefs[0] = 1.0;
bCoefs[1] = D1 / D0;
bCoefs[2] = 0.0;
gain *= (N0 / D0);
start = 1;
}
// Handle quadratic factors. j indexes quad factors
// since each quad factor has three coefs.
for( i=start; i<numQuads; ++i )
{
j = 3 * i;
N0 = aCoefs[j]*f4 + aCoefs[j+1]*f2 + aCoefs[j+2];
N1 = 2 * ( aCoefs[j+2] - aCoefs[j]*f4 );
N2 = aCoefs[j]*f4 - aCoefs[j+1]*f2 + aCoefs[j+2];
D0 = bCoefs[j]*f4 + bCoefs[j+1]*f2 + bCoefs[j+2];
D1 = 2 * ( bCoefs[j+2] - bCoefs[j]*f4 );
D2 = bCoefs[j]*f4 - bCoefs[j+1]*f2 + bCoefs[j+2];
aCoefs[j] = 1.0;
aCoefs[j+1] = N1 / N0;
aCoefs[j+2] = N2 / N0;
bCoefs[j] = 1.0;
bCoefs[j+1] = D1 / D0;
bCoefs[j+2] = D2 / D0;
gain *= (N0 / D0);
}
}
/**
* Unwarps analog frequency for IIR filtrs.
*/
void IIR::freqUnwrap()
{
if( filtType=="lowpass" || filtType=="highpass" )
{
wpass1 = 2*fsamp * atan( wpass1 / (2*fsamp) );
wstop1 = 2*fsamp * atan( wstop1 / (2*fsamp) );
}
else
{
wpass1 = 2*fsamp * atan( wpass1 / (2*fsamp) );
wpass2 = 2*fsamp * atan( wpass2 / (2*fsamp) );
wstop1 = 2*fsamp * atan( wstop1 / (2*fsamp) );
wstop2 = 2*fsamp * atan( wstop2 / (2*fsamp) );
}
}
/**
* Calculates normal Butterworth filter coefficients.
*/
void IIR::butterCoefs()
{
int m, a, b;
double R, epsilon, theta, sigma, omega;
// Make calculations of necessary constants.
epsilon = sqrt( pow(10.0,-0.1*apass1) - 1.0 );
R = pow( epsilon, -1.0/order );
// Initialize gain to 1.0, start indices at 0.
gain = 1.0;
a = 0;
b = 0;
// Handle odd order if necessary.
if( order % 2 )
{
aCoefs[a++] = 0.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = R;
bCoefs[b++] = 0.0;
bCoefs[b++] = 1.0;
bCoefs[b++] = R;
}
// Handle all quadratic terms.
for( m=0; m<order/2; ++m )
{
// Calculate angle first, then real and imag position.
theta = PI*( 2*m + order+1 ) / ( 2*order );
sigma = R * cos(theta);
omega = R * sin(theta);
// Set the quadratic coefficients.
aCoefs[a++] = 0.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = sigma*sigma + omega*omega;
bCoefs[b++] = 1.0;
bCoefs[b++] = -2*sigma;
bCoefs[b++] = sigma*sigma + omega*omega;
}
}
/**
* Calculaets normal Chebyshev coefficients.
*/
void IIR::chebyCoefs()
{
int m, a, b;
double D, epsilon, phi, sigma, omega;
// Make calculations of necessary constants.
epsilon = sqrt( pow(10.0,-0.1*apass1) - 1.0 );
D = asinh(1/epsilon) / order;
// Handle odd order if necessary.
a = 0;
b = 0;
if( order % 2 )
{
gain = 1.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = sinh(D);
bCoefs[b++] = 0.0;
bCoefs[b++] = 1.0;
bCoefs[b++] = sinh(D);
}
else
gain = pow( 10.0, 0.05*apass1 );
// Handle all quadratic terms.
for( m=0; m<order/2; m++ )
{
// Calculate angle first, then real and imag position.
phi = PI * ( 2*m+1 ) / ( 2*order );
sigma = -1 * sinh(D) * sin(phi);
omega = cosh(D) * cos(phi);
// Set the quadratic coefs.
aCoefs[a++] = 0.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = sigma*sigma + omega*omega;
bCoefs[b++] = 1.0;
bCoefs[b++] = -2 * sigma;
bCoefs[b++] = sigma*sigma + omega*omega;
}
}
/**
* Calculates normal inverse Chebyshev coefficients.
*/
void IIR::invChebyCoefs()
{
int m, a, b;
double D, epsilon, mag2, phi, zero, sigma, omega;
// Make calculations of necessary constants.
epsilon = 1 / sqrt( pow(10.0,-0.1*astop1) - 1.0 );
D = asinh( 1/epsilon ) / order;
// Initialize gain to 1.0, start indices at 0.
gain = 1.0;
a = 0;
b = 0;
// Handle odd order if necessary.
if( order % 2 )
{
aCoefs[a++] = 0.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = 1 / sinh(D);
bCoefs[b++] = 0.0;
bCoefs[b++] = 1.0;
bCoefs[b++] = 1 / sinh(D);
}
// Handle all quadratic terms.
for( m=0; m<order/2; m++ )
{
// Calculate angle first, then real and imag position.
phi = PI * ( 2*m+1) / ( 2*order );
sigma = -1 * sinh(D) * sin(phi);
omega = cosh(D) * cos(phi);
// Calculate inverse of pole locations and zero location.
mag2 = omega*omega + sigma*sigma;
omega = -1 * omega / mag2;
sigma = sigma / mag2;
zero = 1. / cos(phi);
// Set the quadratic coefficients.
mag2 = omega*omega + sigma*sigma;
aCoefs[a++] = 1.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = zero * zero;
bCoefs[b++] = 1.0;
bCoefs[b++] = -2 * sigma;
bCoefs[b++] = mag2;
gain *= mag2 / (zero*zero);
}
}
/**
* Calculates normal elliptic coefficients.
*/
void IIR::elliptCoefs()
{
int m, a, b, odd;
double kernel, ratio, epsilon, kn, rt, sigma, omega, zero,
CEIrt, CEIkn, vo, fm,
SP, CP, DP, SN, CN, DN;
// Make calculations of necessary constants.
epsilon = sqrt( pow(10.0,-0.1*apass1) - 1.0 );
// Calculate frequency ratio based on filter selectivity.
if( filtType == "lowpass" )
ratio = wstop1 / wpass1;
else if( filtType == "highpass" )
ratio = wpass1 / wstop1;
else if( filtType == "bandpass" )
ratio = (wstop2-wstop1) / (wpass2-wpass1);
else
ratio = (wpass2-wpass1) / (wstop2-wstop1);
// Calculate values used in future calculations.
kernel = ( pow(10,-0.1*astop1) - 1 ) / ( pow(10,-0.1*apass1) - 1 );
rt = 1 / ratio;
kn = 1 / sqrt(kernel);
// Calculate ellip integrals, vo and ellipFun.
CEIrt = ellipticIntegral(rt);
CEIkn = ellipticIntegral(kn);
vo = ( CEIrt / (CEIkn*order) ) * arcsc( 1/epsilon, kn );
ellipticFun( vo, sqrt(1-rt*rt), &SP, &CP, &DP );
// Handle odd order if necessary.
a = 0;
b = 0;
odd = order % 2;
if( odd )
{
gain = 1.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = SP * CP / (1-SP*SP);
bCoefs[b++] = 0.0;
bCoefs[b++] = 1.0;
bCoefs[b++] = SP * CP / (1-SP*SP);
}
else
gain = pow( 10.0, 0.05*apass1 );
// Handle all quadratic terms.
for( m=0; m<order/2; ++m )
{
// Make intermediate calculations.
fm = CEIrt * ( 2*m+1+odd ) / order;
ellipticFun( fm, rt, &SN, &CN, &DN );
// Calculate real and imag coordinates of poles.
sigma = -1*CN*DN*SP*CP / ( 1 - DN*DN*SP*SP );
omega = SN*DP / ( 1 - DN*DN*SP*SP );
zero = 1 / (rt*SN);
// Set the quadratic coefficients.
aCoefs[a++] = 1.0;
aCoefs[a++] = 0.0;
aCoefs[a++] = zero*zero;
bCoefs[b++] = 1.0;
bCoefs[b++] = -2*sigma;
bCoefs[b++] = sigma*sigma + omega*omega;
gain *= (sigma*sigma + omega*omega) / (zero*zero);
}
}
/**
* Converts normal lowpass coefs to unnormal lowpass
* coefficients at a specific frequency.
*/
void IIR::unnormLP( double freq )
{
int qd, cf, qdStart;
// Handle first order, if odd; set qdStart.
if( order % 2 )
{
aCoefs[2] *= freq;
bCoefs[2] *= freq;
qdStart = 1;
}
else
qdStart = 0;
// Handle quadratic factors, qd indexes through
// quadratic factors, cf converts to coef number.
for( qd=qdStart; qd<(order+1)/2; ++qd )
{
cf = qd * 3;
aCoefs[cf+1] *= freq;
aCoefs[cf+2] *= (freq*freq);
bCoefs[cf+1] *= freq;
bCoefs[cf+2] *= (freq*freq);
}
}
/**
* Converts normal lowpass coefs to unnormal highpass
* coefficients at a specific frequency.
*/
void IIR::unnormHP( double freq )
{
int qd, cf, qdStart;
// Handle first order, if odd; set qdStart
if( order % 2 )
{
gain *= ( aCoefs[2] / bCoefs[2] );
aCoefs[2] = freq*aCoefs[1] / aCoefs[2];
aCoefs[1] = 1.0;
bCoefs[2] = freq / bCoefs[2];
qdStart = 1;
}
else
qdStart = 0;
// Handle quadratic factors, qd indexes through
// quadratic factors, cf converts to coef number.
for( qd=qdStart; qd<(order+1)/2; ++qd )
{
cf = qd * 3;
gain *= ( aCoefs[cf+2] / bCoefs[cf+2] );
aCoefs[cf+1] *= ( freq / aCoefs[cf+2] );
aCoefs[cf+2] = freq * freq *aCoefs[cf]
/ aCoefs[cf+2];
aCoefs[cf] = 1.0;
bCoefs[cf+1] *= ( freq / bCoefs[cf+2] );
bCoefs[cf+2] = freq * freq * bCoefs[cf]
/ bCoefs[cf+2];
bCoefs[cf] = 1.0;
}
}
/**
* Converts normal lowpass coefficients to unnormal
* bandpass coefficients at a specific frequency.
*/
void IIR::unnormBP( double BW, double Wo )
{
int qd, ocf, ncf, qdStart,
numCoefs, orgQuads, orgOrder;
complex<double> A, B, C, D, E;
Vector<double> orgNum,
orgDen,
newNum,
newDen;
// Store orig number of quadratics and order,
// new order will be twice original.
orgOrder = order;
orgQuads = (orgOrder + 1)/2;
order = orgOrder * 2;
orgNum = aCoefs;
orgDen = bCoefs;
numCoefs = 3 * orgOrder;
newNum.resize(numCoefs);
newDen.resize(numCoefs);
// If orgOrder odd, convert first order factor to
// quadratic, qdStart indic start pt for loop.
if( orgOrder % 2 )
{
newNum[0] = orgNum[1];
newNum[1] = BW * orgNum[2];
newNum[2] = orgNum[1] * Wo * Wo;
newDen[0] = orgDen[1];
newDen[1] = BW * orgDen[2];
newDen[2] = orgDen[1] * Wo * Wo;
qdStart = 1;
}
else
qdStart = 0;
// Each orig quad term will be converted to two new
// quadratics via complex quadratic factoring.
for( qd=qdStart; qd<orgQuads; ++qd )
{
ocf = qd * 3;
ncf = qd * 6 - qdStart * 3;
// For numers which DON'T have s^2 or s terms.
if( orgNum[ocf] == 0.0 )
{
newNum[ncf] = 0.0;
newNum[ncf+1] = sqrt(orgNum[ocf+2]) * BW;
newNum[ncf+2] = 0.0;
newNum[ncf+3] = 0.0;
newNum[ncf+4] = sqrt(orgNum[ocf+2]) * BW;
newNum[ncf+5] = 0.0;
}
// For numers which DO have s^2 and s terms.
else
{
A = complex<double>( orgNum[ocf], 0.0 );
B = complex<double>( orgNum[ocf+1], 0.0 );
C = complex<double>( orgNum[ocf+2], 0.0 );
quadradicRoot( A, B, C, D, E );
A = complex<double>( 1.0, 0.0 );
B = - D * complex<double>( BW, 0.0 );
C = complex<double>( Wo*Wo, 0.0 );
quadradicRoot( A, B, C, D, E );
newNum[ncf] = 1.0;
newNum[ncf+1] = -2.0 * real(D);
newNum[ncf+2] = norm(D);
newNum[ncf+3] = 1.0;
newNum[ncf+4] = -2.0 * real(E);
newNum[ncf+5] = norm(E);
}
A = complex<double>( orgDen[ocf], 0.0 );
B = complex<double>( orgDen[ocf+1], 0.0 );
C = complex<double>( orgDen[ocf+2], 0.0 );
quadradicRoot( A, B, C, D, E );
// Make required substitutions, factor again.
A = complex<double>( 1.0, 0.0 );
B = - D * complex<double>( BW, 0.0 );
C = complex<double>( Wo*Wo, 0.0 );
quadradicRoot( A, B, C, D, E );
// Make required substitutions, factor again.
newDen[ncf] = 1.0;
newDen[ncf+1] = -2.0 * real(D);
newDen[ncf+2] = norm(D);
newDen[ncf+3] = 1.0;
newDen[ncf+4] = -2.0 * real(E);
newDen[ncf+5] = norm(E);
}
aCoefs = newNum;
bCoefs = newDen;
}
/**
* Converts normal lowpass coefficients to unnormal
* bandstop coefficients at a specific frequency.
*/
void IIR::unnormBS( double BW, double Wo )
{
int qd, ocf, ncf, qdStart,
numCoefs, orgQuads, orgOrder;
complex<double> A, B, C, D, E;
Vector<double> orgNum,
orgDen,
newNum,
newDen;
// Store orig number of quadratics and order,
// new order will be twice original.
orgOrder = order;
orgQuads = (orgOrder + 1)/2;
order = orgOrder * 2;
orgNum = aCoefs;
orgDen = bCoefs;
numCoefs = 3 * orgOrder;
newNum.resize(numCoefs);
newDen.resize(numCoefs);
// If orgOrder odd, convert first order factor to
// quadratic, qdStart indic start pt for loop.
if( orgOrder % 2 )
{
gain *= (orgNum[2] / orgDen[2]);
newNum[0] = 1.0;
newNum[1] = BW * orgNum[1] / orgNum[2];
newNum[2] = Wo * Wo;
newDen[0] = 1.0;
newDen[1] = BW * orgDen[1] / orgDen[2];
newDen[2] = Wo * Wo;
qdStart = 1;
}
else
qdStart = 0;
// Each orig quad term will be converted to two new
// quadratics via complex quadratic factoring.
for( qd=qdStart; qd<orgQuads; qd++ )
{
ocf = qd * 3;
ncf = qd * 6 - qdStart * 3;
gain *= (orgNum[ocf+2] / orgDen[ocf+2]);
// For numers which DON'T have s^2 or s terms.
if(orgNum[ocf] == 0)
{
newNum[ncf] = 1.0;
newNum[ncf+1] = 0.0;
newNum[ncf+2] = Wo * Wo;
newNum[ncf+3] = 1.0;
newNum[ncf+4] = 0.0;
newNum[ncf+5] = Wo * Wo;
}
// For numers which DO have s^2 and s terms.
else
{
A = complex<double>( orgNum[ocf], 0.0 );
B = complex<double>( orgNum[ocf+1], 0.0 );
C = complex<double>( orgNum[ocf+2], 0.0 );
quadradicRoot( A, B, C, D, E );
A = complex<double>( 1.0, 0.0 );
B = - 1.0/D * complex<double>( BW, 0.0 );
C = complex<double>( Wo*Wo, 0.0 );
quadradicRoot( A, B, C, D, E );
newNum[ncf] = 1.0;
newNum[ncf+1] = -2.0 * real(D);
newNum[ncf+2] = norm(D);
newNum[ncf+3] = 1.0;
newNum[ncf+4] = -2.0 * real(E);
newNum[ncf+5] = norm(E);
}
A = complex<double>( orgDen[ocf], 0.0 );
B = complex<double>( orgDen[ocf+1], 0.0 );
C = complex<double>( orgDen[ocf+2], 0.0 );
quadradicRoot( A, B, C, D, E );
// Make required substitutions, factor again.
A = complex<double>( 1.0, 0.0 );
B = - 1.0/D * complex<double>( BW, 0.0 );
C = complex<double>( Wo*Wo, 0.0 );
quadradicRoot( A, B, C, D, E );
// Make required substitutions, factor again.
newDen[ncf] = 1.0;
newDen[ncf+1] = -2.0 * real(D);
newDen[ncf+2] = norm(D);
newDen[ncf+3] = 1.0;
newDen[ncf+4] = -2.0 * real(E);
newDen[ncf+5] = norm(E);
}
aCoefs = newNum;
bCoefs = newDen;
}
/**
* Calculates response for IIR filters.
*/
double IIR::frqeResp( double freq )
{
int c, q;
double omega, omega2,
rea, img, tmp,
mag;
// Initialize magnatude.
mag = gain;
// Pre-calculate omega and omega squared.
omega = TWOPI * freq / fsamp;
omega2 = 2 * omega;
// Loop through coefficients for each quadratic.
for( q=0; q<(order+1)/2; q++ )
{
c = 3*q;
rea = aCoefs[c] + aCoefs[c+1]*cos(omega) + aCoefs[c+2]*cos(omega2);
img = -aCoefs[c+1]*sin(omega) - aCoefs[c+2]*sin(omega2);
tmp = sqrt( rea*rea + img*img );
mag *= tmp;
rea = bCoefs[c] + bCoefs[c+1]*cos(omega) + bCoefs[c+2]*cos(omega2);
img = -bCoefs[c+1]*sin(omega) - bCoefs[c+2]*sin(omega2);
tmp = sqrt(rea*rea + img*img);
mag /= tmp;
}
// Convert magnitude response to dB.
if( mag < EPS )
mag = EPS;
mag = 20 * log10(mag);
return mag;
}
/**
* Calculate the response at edge freqency.
*/
void IIR::calcGain()
{
// Determine the edge frequency.
if( filtType=="lowpass" || filtType=="highpass" )
{
double f1 = wpass1 / TWOPI,
f2 = wstop1 / TWOPI;
edgeGain.resize(2);
edgeGain(1) = frqeResp( f1 );
edgeGain(2) = frqeResp( f2 );
}
else
{
double f1 = wpass1 / TWOPI,
f2 = wpass2 / TWOPI,
f3 = wstop1 / TWOPI,
f4 = wstop2 / TWOPI;
edgeGain.resize(4);
edgeGain(1) = frqeResp( f1 );
edgeGain(2) = frqeResp( f2 );
edgeGain(3) = frqeResp( f3 );
edgeGain(4) = frqeResp( f4 );
}
}测试代码：
/*****************************************************************************
* iir_test.cpp
*
* IIR class testing.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iir.h>
using namespace std;
using namespace splab;
int main()
{
// string fType = "lowpass";
// double fs = 1000,
// fpass = 200,
// apass = -3,
// fstop = 300,
// astop = -20;
// string method = "Butterworth";
// string method = "Chebyshev";
// string method = "InvChebyshev";
// string method = "Elliptic";
// IIR iir( fType, method );
// iir.setParams( fs, fpass, apass, fstop, astop );
// string fType = "highpass";
// double fs = 1000,
// fstop = 200,
// astop = -20,
// fpass = 300,
// apass = -3;
// string method = "Butterworth";
// string method = "Chebyshev";
// string method = "InvChebyshev";
// string method = "Elliptic";
// IIR iir( fType, method );
// iir.setParams( fs, fstop, astop, fpass, apass );
// string fType = "bandpass";
// double fs = 1000,
// fstop1 = 100,
// astop1 = -20,
// fpass1 = 200,
// fpass2 = 300,
// apass1 = -3,
// fstop2 = 400,
// astop2 = -20;
// string method = "Butterworth";
// string method = "Chebyshev";
// string method = "InvChebyshev";
// string method = "Elliptic";
// IIR iir( fType, method );
// iir.setParams( fs, fstop1, astop1, fpass1, fpass2, apass1, fstop2, astop2 );
string fType = "bandstop";
double fs = 1000,
fpass1 = 100,
apass1 = -3,
fstop1 = 200,
fstop2 = 300,
astop1 = -20,
fpass2 = 400,
apass2 = -3;
// string method = "Butterworth";
// string method = "Chebyshev";
// string method = "InvChebyshev";
string method = "Elliptic";
IIR iir( fType, method );
iir.setParams( fs, fpass1, apass1, fstop1, fstop2, astop1, fpass2, apass2 );
iir.design();
iir.dispInfo();
cout << endl;
return 0;
}运行结果：
Filter selectivity : bandstop
Approximation method : Elliptic
Filter order : 4
Overall gain : 0.153272
Sampling Frequency (Hz) : 1000
Lower passband frequency (Hz) : 100
Lower passband gain (dB) : -3
Lower stopband frequency (Hz) : 200
Upper stopband frequency (Hz) : 300
Stopband gain (dB) : -20
Upper passband frequency (Hz) : 400
Upper passband gain (dB) : -3
Numerator Coefficients Denominator Coefficients
N [ 1 z^-1 z^-2 ] [ 1 z^-1 z^-2 ]
=== ================================= =================================
1 1.0 -0.45087426 1.00000000 1.0 -1.44482415 0.70572930
2 1.0 0.45087426 1.00000000 1.0 1.44482415 0.70572930
==================== Edge Frequency Response ====================
Mag(fp1) = -3.00000000(dB) Mag(fp2) = -3.00000000(dB)
Mag(fs1) = -36.87418446(dB) Mag(fs2) = -36.87418446(dB)
Process returned 0 (0x0) execution time : 0.078 s
Press any key to continue.
扫一扫
727
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
