头文件：
/*
* 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
*/
/*****************************************************************************
* dfd.h
*
* Digital Filter Design.
*
* This is a base class for FIR and IIR filter design. It contains the design
* specifies and some universal parameters of filter.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef DFD_H
#define DFD_H
#include <string>
#include <vector.h>
namespace splab
{
const double FREQ_MIN = 1.0E-3; // min freqency
const double FREQ_MAX = 1.0E12; // max freqency
const double GAIN_PASS = -1.0E-2; // max passband gain
const double GAIN_TRAN = -3.0103; // min pb, max sb gain
const double GAIN_STOP = -1.0E02; // min stopband gain
class DFD
{
public:
DFD( const string &select );
virtual ~DFD();
void setParams( double fs, double f1, double a1,
double f2, double a2 );
void setParams( double fs, double f1, double a1, double f2,
double f3, double a2, double f4, double a3 );
virtual void design() = 0;
virtual void dispInfo() const = 0;
protected:
double fsamp; // sampling frequency
double wpass1, wpass2; // passband edge frequency
double wstop1, wstop2; // stopband edge frequency
double apass1, apass2; // passband gain
double astop1, astop2; // stopband gain
string filtType; // filter type
int order; // length of filter
inline double getValue( double x, double min, double max );
};
// class DFD
#include <dfd-impl.h>
}
// namespace splab
#endif
// DFD_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
*/
/*****************************************************************************
* fir.h
*
* Finite Impulse Response Digital Filter.
*
* This class is designed for designing the FIR filter using window function
* 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 windows are:
* Rectangle Bartlett Blackman Hanning Hamming
* Gauss Kaiser
*
* The length of filter( filter's order plus one ) are multiples of 4, which
* is the least number(L=4n) satisfying the design specifies.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef FIR_H
#define FIR_H
#include <iomanip>
#include <window.h>
#include <dfd.h>
namespace splab
{
using std::setw;
using std::ios;
using std::setiosflags;
using std::setprecision;
class FIR : public DFD
{
public:
FIR( const string &select, const string &win );
FIR( const string &select, const string &win, double a );
~FIR();
void design();
void dispInfo() const;
Vector<double> getCoefs() const;
private:
void orderEst();
void idealCoef();
void calcCoef();
double frqeResp( double freq );
void calcGain();
bool isSatisfy();
string windType; // window type
Vector<double> wind, // window function
coefs, // coefficients
edgeGain; // gain at edge frquency
double alpha; // window parameter
};
// class FIR
#include <fir-impl.h>
}
// namespace splab
#endif
// FIR_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
*/
/*****************************************************************************
* dfd-impl.h
*
* Implementation for DFD class.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
/**
* constructors and destructor
*/
DFD::DFD( const string &select ) : filtType(select)
{
}
DFD::~DFD()
{
}
/**
* Get a number which falls between min and max.
*/
double DFD::getValue( double x, double min, double max )
{
if( (x > min) && (x < max) )
return x;
else
{
cerr << "The parameter is out of range!" << endl;
return 0;
}
}
/**
* Set filter's parameters
*/
void DFD::setParams( double fs, double f1, double a1, double f2, double a2 )
{
fsamp = getValue( fs, FREQ_MIN, FREQ_MAX );
double maxFreq = fsamp / 2;
if( filtType == "lowpass" )
{
wpass1 = getValue( f1, FREQ_MIN, maxFreq );
apass1 = getValue( a1, GAIN_TRAN, GAIN_PASS );
wstop1 = getValue( f2, wpass1, maxFreq );
astop1 = getValue( a2, GAIN_STOP, GAIN_TRAN );
}
else if( filtType == "highpass" )
{
wstop1 = getValue( f1, FREQ_MIN, maxFreq );
astop1 = getValue( a1, GAIN_STOP, GAIN_TRAN );
wpass1 = getValue( f2, wstop1, maxFreq );
apass1 = getValue( a2, GAIN_TRAN, GAIN_PASS );
}
else
cerr << "Parameters setting has failed!" << endl;
// Convert all edge frequencies to rad/sec
wpass1 *= TWOPI;
wstop1 *= TWOPI;
}
void DFD::setParams( double fs, double f1, double a1, double f2,
double f3, double a2, double f4, double a3 )
{
fsamp = getValue( fs, FREQ_MIN, FREQ_MAX );
double maxFreq = fsamp / 2;
if( filtType == "bandpass" )
{
wstop1 = getValue( f1, FREQ_MIN, maxFreq );
astop1 = getValue( a1, GAIN_STOP, GAIN_TRAN );
wpass1 = getValue( f2, wstop1, maxFreq );
wpass2 = getValue( f3, wpass1, maxFreq );
apass1 = getValue( a2, GAIN_TRAN, GAIN_PASS );
wstop2 = getValue( f4, wpass2, maxFreq );
astop2 = getValue( a3, GAIN_STOP, GAIN_TRAN );
apass2 = apass1;
if( astop1 < astop2 )
astop2 = astop1;
else
astop1 = astop2;
}
else if( filtType == "bandstop" )
{
wpass1 = getValue( f1, FREQ_MIN, maxFreq );
apass1 = getValue( a1, GAIN_TRAN, GAIN_PASS );
wstop1 = getValue( f2, wpass1, maxFreq );
wstop2 = getValue( f3, wstop1, maxFreq );
astop1 = getValue( a2, GAIN_STOP, GAIN_TRAN );
wpass2 = getValue( f4, wstop2, maxFreq );
apass2 = getValue( a3, GAIN_TRAN, GAIN_PASS );
astop2 = astop1;
if( apass1 < apass2 )
apass1 = apass2;
else
apass2 = apass1;
}
else
cerr << "Parameters setting has failed!" << endl;
// Convert all edge frequencies to rad/sec
wpass1 *= TWOPI;
wpass2 *= TWOPI;
wstop1 *= TWOPI;
wstop2 *= TWOPI;
}/*
* 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
*/
/*****************************************************************************
* fir-impl.h
*
* Implementation for FIR class.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
/**
* constructors and destructor
*/
FIR::FIR( const string &select, const string &win )
: DFD( select ), windType(win)
{
bool cond = ( filtType=="lowpass" || filtType=="highpass" ||
filtType=="bandpass" || filtType=="bandstop" );
assert(cond);
cond = ( windType=="Rectangle" || windType=="Bartlett" ||
windType=="Blackman" || windType=="Hanning" ||
windType=="Hamming" || windType=="Gauss" );
assert(cond);
}
FIR::FIR( const string &select, const string &win, double a )
: DFD( select ), windType(win), alpha(a)
{
bool cond = ( filtType=="lowpass" || filtType=="highpass" ||
filtType=="bandpass" || filtType=="bandstop" );
assert(cond);
cond = ( windType=="Kaiser" || windType=="Gauss" );
assert(cond);
}
FIR::~FIR()
{
}
/**
* Design digital FIR filter.
*/
void FIR::design()
{
// Estimate the length of filter.
orderEst();
calcCoef();
calcGain();
while( !isSatisfy() )
{
order += 4;
calcCoef();
calcGain();
}
}
/**
* Displays the design result of filter.
*/
void FIR::dispInfo() const
{
int i, j, k;
cout << endl;
cout << "\t\t Filter selectivity : " << filtType << endl;
cout << "\t\t Window type : " << windType << 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;
cout << "\t\t\t\t Filter Coefficients" << endl << endl;
cout << " N [ N + 0 N + 1" ;
cout << " N + 2 N + 3 ]" << endl;
cout << " === ============================" ;
cout << "========================================" ;
for( i=0; i<order/4; ++i )
{
j = i * 4;
cout << endl << setw(4) << j << " ";
cout << setiosflags(ios::scientific) << setprecision(8);
for( k=0; k<4; ++k )
cout << setw(16) << coefs[j+k] << " " ;
}
cout << endl << 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> FIR::getCoefs() const
{
return coefs;
}
/**
* Estimates the length of FIR filter.
*/
void FIR::orderEst()
{
double deltaFreq, // delta frequency
lowerDF, // delta freq of lower band
upperDF, // delta freq of upper band
errSB, // stopband error
errPB, // passband error
errMin, // minimum of sb/pb errors
errDB; // minimum error in dB
// Determine frequency delta
if( filtType=="lowpass" || filtType=="highpass" )
deltaFreq = abs(wstop1-wpass1) / fsamp;
else
{
lowerDF = abs(wstop1-wpass1) / fsamp;
upperDF = abs(wstop2-wpass2) / fsamp;
if( lowerDF > upperDF )
deltaFreq = upperDF;
else
deltaFreq = lowerDF;
}
// Determine stopband and passband errors
errPB = 1 - pow( 10, 0.05*apass1 );
errSB = pow( 10, 0.05*astop1 );
if( errSB < errPB )
errMin = errSB;
else
errMin = errPB;
errDB = -20 * log10(errMin);
// Store filter length in pFilt and return beta.
if( errDB > 21 )
order = int( ceil( 1 + (errDB-7.95) / (2.285*deltaFreq) ) );
else
order = int( ceil( 1 + (5.794/deltaFreq) ) );
order += 4 - order%4;
}
/**
* Calculates ideal FIR coefficients.
*/
void FIR::idealCoef()
{
int i;
double t,
tau,
Wc1,
Wc2;
// Calculate tau as non-integer if order is even.
tau = ( order - 1 ) / 2.0;
// Adjust cutoff frequency to midway point between stop
// and pass edge frequency and convert to digital frequency.
Wc1 = (wstop1 + wpass1) / (2*fsamp);
Wc2 = (wstop2 + wpass2) / (2*fsamp);
// Calc coefs based on selectivity of filter.
if( filtType == "lowpass" )
for( i=0; i<order; ++i )
{
if( i == tau )
coefs[i] = Wc1/PI;
else
{
t = i - tau;
coefs[i] = sin(Wc1*t) / (PI*t);
}
}
else if( filtType == "highpass" )
for( i=0; i<order; ++i )
{
if( i == tau )
coefs[i] = (PI-Wc1) / PI;
else
{
t = i - tau;
coefs[i] = ( sin(PI*t) - sin(Wc1*t) ) / (PI*t);
}
}
else if( filtType == "bandpass" )
for( i=0; i<order; ++i )
{
if( i == tau )
coefs[i] = ( Wc2-Wc1 ) / PI;
else
{
t = i - tau;
coefs[i] = ( sin(Wc2*t) - sin(Wc1*t) ) / (PI*t);
}
}
else
for( i=0; i<order; ++i )
{
if( i == tau )
coefs[i] = ( PI+Wc1-Wc2 ) / PI;
else
{
t = i - tau;
coefs[i] = ( sin(PI*t)-sin(Wc2*t)+sin(Wc1*t) ) / (PI*t);
}
}
}
/**
* Calculates the digital FIR coefficients.
*/
void FIR::calcCoef()
{
coefs.resize(order);
wind.resize(order);
// Calculate the ideal FIR coefficients.
idealCoef();
// Get window function.
if( windType=="Kaiser" || windType=="Gauss" )
wind = window( windType, order, alpha, 1.0 );
else
wind = window( windType, order, 1.0 );
// Multiply window and ideal coefficients.
coefs *= wind;
}
/**
* Calculate the response at edge freqency.
*/
double FIR::frqeResp( double freq )
{
double mag = 1.0,
rea = 0.0,
img = 0.0,
omega = TWOPI*freq / fsamp;
// Loop through all the coefficients.
for( int i=0; i<order; ++i )
{
double domega = i * omega;
rea += coefs[i] * cos(domega);
img += coefs[i] * sin(domega);
}
// Calculate final result and convert to degrees.
mag = sqrt( rea*rea + img*img );
if( mag < EPS )
mag = EPS;
mag = 20 * log10( mag );
return mag;
}
/**
* Calculate the response at edge freqency.
*/
void FIR::calcGain()
{
// Determine the edge freqency.
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 );
}
}
/**
* Design digital FIR filter.
*/
bool FIR::isSatisfy()
{
if( edgeGain.size() == 2 )
{
if( edgeGain(1)<apass1 || edgeGain(2)>astop1 )
return false;
}
else
{
if( edgeGain(1)<apass1 || edgeGain(2)<apass2 ||
edgeGain(3)>astop1 || edgeGain(4)>astop2 )
return false;
}
return true;
}测试代码：
/*****************************************************************************
* fir_test.cpp
*
* FIR class testing.
*
* Zhang Ming, 2010-03, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <fir.h>
using namespace std;
using namespace splab;
int main()
{
string wType = "Hamming";
// string fType = "lowpass";
// double fs = 1000,
// fpass = 200,
// apass = -3,
// fstop = 300,
// astop = -20;
// FIR fir( fType, wType );
// fir.setParams( fs, fpass, apass, fstop, astop );
// string fType = "highpass";
// double fs = 1000,
// fstop = 200,
// astop = -20,
// fpass = 300,
// apass = -3;
// FIR fir( fType, wType );
// fir.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;
// FIR fir( fType, wType );
// fir.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;
FIR fir( fType, wType );
fir.setParams( fs, fpass1, apass1, fstop1, fstop2, astop1, fpass2, apass2 );
fir.design();
fir.dispInfo();
cout << endl;
return 0;
}Filter selectivity : bandstop
Window type : Hamming
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
Filter Coefficients
N [ N + 0 N + 1 N + 2 N + 3 ]
=== ====================================================================
0 -3.85192764e-003 8.42472981e-003 3.11153775e-003 -2.03549729e-003
4 -3.55185234e-002 2.30171390e-002 -1.41331145e-001 2.61755439e-001
8 2.88881794e-002 3.56158158e-001 3.56158158e-001 2.88881794e-002
12 2.61755439e-001 -1.41331145e-001 2.30171390e-002 -3.55185234e-002
16 -2.03549729e-003 3.11153775e-003 8.42472981e-003 -3.85192764e-003
==================== Edge Frequency Response ====================
Mag(fp1) = -0.85449456(dB) Mag(fp2) = -0.80635855(dB)
Mag(fs1) = -21.347821(dB) Mag(fs2) = -21.392825(dB)
Process returned 0 (0x0) execution time : 0.078 s
Press any key to continue.
939
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
