/*
FIR 滤波器设计之切比雪夫(Chebyshev)逼近法
void defir3(int nfilt, int nbands, float edge[], float fx[], float wtx[], float h[])
输入：
nfilt FIR滤波器长度(单位抽样响应的长度),< LenFilt
nbands 频带数(通带和阻带的总数目),<= MaxBands
fx 函数数组(fx[nbands]),各个频带的理想频率响应幅值(如通带设为1, 阻带设为0)
wtx 权函数数组(wtx[nbands]),每个频带的频率抽样权重数组(控制通带和阻带的衰减)
edge 频带边沿数组(通带阻带的带边归一化频率),频率下边界与上边界(edge[2*nbands])
输出：
h 单位抽样响应 h[nfilt]
代码来源及参考文献:
[1] 胡广书 编著《数字信号处理——理论、算法与实现》1997.8 清华大学出版社 P.250,471
加权切比雪夫最佳一致逼近FIR设计法
[2] H.McCLELLAN, A computer program for designing optimum FIR linear phase digital filters,
IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS, VOL.AU-21, NO.6, DECEMBER 1973
[3] John G. Proakis 等, 数字信号处理：理论、算法与应用 (第三版),电子工业出版社, 2004.5, P.522~536
说明
1) 此C代码最早来自对[1]中Fortran程序的改编(并简化了部分语句行和流程), [1]已省略[2]中差分和Hilbert变换内容.
后来参考[2]流程图及代码, 加入了[2]中的注释(每段前与代码行齐平), 并附加个人理解或疑问(行末及缩进注释), 及些许修改.
2) 对于研究该算法, [2]是必备的, 因为[1]中有些许错误且不详细([1]1997.8第一版P.253),
但[2]中也可能有错误, 如函数d()的算法导致迭代结果不正确; P.7/21(P.512) Fig.5中第1个"DEV<=~ERR"分支Yes/No反了.
3) 本程序已调通, 但可读性仍很差, 尤其Remez算法中搜索局部极值点的过程,
[2]的流程居然超乎想象的繁琐, 似乎不应该, 有待研究.
4) 本文档主要供对原始算法研究, 注释繁杂, 也有个人的疑问. 欢迎各位同仁指教或交流心得.
——2009.9.5~9 liu-qg@tom.com
*/
#include <math.h>
#include <stdio.h>
bool debugView =true;
int chView =0;
//#define I_have_VDAA 1
#ifdef I_have_VDAA
//应用外部独立程序 VDAA.exe 显示曲线细节及其动态变化:
// 将迭代过程中各种中间结果发送到VDAA
#include "..\..\VDAA\IPC\VDAA.h"
VDAA vdaa;
#define vdaa_putFloat vdaa.putFloat
#define vdaa_put1D vdaa.put1D
#define vdaa_ready vdaa.ready
#define vdaa_setAutoRefresh vdaa.setAutoRefresh
//疏集上点值投射到密集
template <class T>
void put2Vdaa(T xm[], int ik[], int M, int N,
int chiBase0, const char *ChName=0, const char *Unit=0)
{
float *tmp =new float [N];
int i=0;
while (i<N) tmp[i++] =0;
i =0;
while (i<M){
tmp[ik[i]-1] =xm[i]; //ik[i]表示的下标从1开始, tmp[]下标从0开始
i++;
}
vdaa_putFloat(tmp, N, chiBase0, 0, ChName, Unit);
delete tmp;
}
#else //I don't have VDAA.exe
//将涉及vdaa的代码全部注释, 不影响运算结果.
#define put2Vdaa
#define vdaa_putFloat
#define vdaa_put1D
#define vdaa_ready()
#define vdaa_setAutoRefresh()
#endif
bool Equal(float a, float b){
const float eps =1.e-12;
float d =a-b;
return (-eps < d && d < eps);
}
const int LGrid =16; //误差函数E(w)插值的网格密度——[2] P.529
const int LenFilt =128; //滤波器长度最大限制值
const int N =LGrid*LenFilt+5;
const int MaxBands =10;
float pi2, dev;
float grid[N], des[N], wt[N]; //在频带密格点处的频率、理想响应、拟合偏差的权重
float Hr[N]; //增加此拟合函数供调试观察
float Error[N]; //增加此目标偏差函数, 便于简化语句
int nfcns; //逼近函数的数量, 与滤波器类型和滤波器长度有关——各引用处的准确含义有待仔细研究...
int ngrid; //对频率点加密后的总格点数
//打算按fortran习惯取用[1~66] fortran规则: [i:j]==[i~j]即[下界:上界],[n]==[1~n]
float x[66+1], y[66+1], alpha[66+1];
double ad[66+1];
int iext[66+1]; //极值频率格点索引，如弟j个极值点频率为grid[iext[j]]
//注意, 以下代码来自 Fortran, 其数组元素起始下标习惯为 1, 改编为C语言后仍保持原样.
void defir3(int nfilt, int nbands, float edge_[], float fx_[], float wtx_[], float h_[])
{
//输入C语法习惯的数组, 以Fortran语法习惯来引用. 二者起始元素在共同的地址:
//
// Fortran--> edge[1] == edge_[0] <--C
//
//故取C语法习惯的数组[-1]元素地址, 视为Fortran习惯的数组,
// 则引用数组元素的Fortran语句可保持原样.
float *edge =edge_-1,
*fx =fx_ -1,
*wtx =wtx_ -1,
*h =h_ -1;
//于是得到Fortran习惯数组 float edge[1:2*nbands], float fx[1:nbands], float wtx[1:nbands], float h[1:nfilt])
pi2 =8.*atan(1.);
float pi =pi2/2.;
float nfmax =LenFilt;
if (MaxBands<nbands) return ;
if (3<=nfilt && nfilt<=nfmax) goto L115;
return ;
L115:
if (nbands<=0) nbands=1;
//nfilt =nfilt +1; [1]中语句, 弃用
//准备工作 [2] Fig.2-----P.4
//逼近函数数量 nfcns =(滤波器长度一半,包括中心点)
int nodd =(nfilt % 2); //nfilt为奇数
nfcns =nfilt/2 +(nfilt % 2);
//[2]中若 NEG=1 表示FIR对称形式为负,如: h[n] =-h[M-n]
//在此假定 NEG=0, 如 h[n] =h[M-n]
//SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID IS
// (FILTER LENGTH + 1)*GRID DENSITY/2
grid[1] =edge[1]; //edge[1],grid[1]是否可以不为频域下界0?
float delf =0.5/(LGrid*nfcns); //密格点间距——只是取大概的值而已, 以便得到适当数量的格点
if (edge[1]<delf) grid[1] =delf; //? grid[首点] != 0, 为何特别剔除频域[0, 0.5]的左边界?
int j =1, l=1, lband=1;
//在频率轴上各个频带按增量 delf=0.5/(lgrid*nfcns) 填充 grid[] (grid[j]=grid[j-1]+delf);
while (true)
{
//L140:
float fup =edge[l+1];
//L145:
//对每个格点计算理想幅频响应函数des[],权重wtx[]
do {
des[j] =fx [lband];
wt [j] =wtx[lband];
grid[j+1] =grid[j]+delf;
j++;
}while (grid[j]<=fup); //取得每个频带上的 des[], wt[], grid[]
//L150:
grid[j-1] =fup; //处理频带端点: 频带边界点强制设为格点上!
//可见格点并非都是严格等间距的,后果是总的个点数 ngrid 事先不能确定!
//des [j-1] =fx [lband]; 上面循环内其实已设置了
//wt [j-1] =wtx[lband];
lband++, l+=2; //跳过了过渡带
if (nbands<lband){
ngrid =j-1; //goto L160
break;
}
grid[j] =edge[l]; //至此可见过渡带内不设格点, 只在指定频带内设格点(包括频带2端点), 除带沿处外格点均匀分布.
//拟合时仅考虑指定频带内的偏差, 过渡带自然成形
}//goto L140;
//L160:
if (nodd==0) //IF(NEG.NE.NODD) GO TO 165
if ((.5-delf)<grid[ngrid]) ngrid =ngrid-1; //?去掉靠近右端0.5的边界点,是何道理. 前面还去掉了左界0
if (debugView)
vdaa_putFloat(grid+1, ngrid, chView, 0, "grid", "");
//L165:
//用一个新的逼近问题等效原逼近问题
if (nodd==1) goto L200;
//change des and wt ... ? 为何 nfilt偶数时 ——参考 [2] P.4, Fig.2
for (j=1; j<=ngrid; j++){
float change =cos(pi*grid[j]); //原文中Fortran函数 DCOS(x) means Double precision Cosine
des[j] /=change; //这里或许解释了前面(nodd==0时)为何去掉0.5附近的最后1格点grid[ngrid]
wt[j] *=change;
} //175
L200:
if (debugView){
vdaa_putFloat(des+1, ngrid, chView++, 0, "des", "");
vdaa_putFloat(wt+1, ngrid, chView, 0, "wt", "");
}
//设置初始极值点频率位置(猜测值)——对格点"等间距"分布, 其索引保存到iext[]
float dg =float(ngrid-1)/float(nfcns); //从ngrid个格点中等间隔取nfcns+1个极值点
for (j=1; j<=nfcns; j++){
iext[j] =(j-1)*dg+1;
}//210
iext[nfcns+1] =ngrid; //将2端点grid[1],grid[ngrid]设为极值点
int nm1 =nfcns -1;
int nz =nfcns +1;
//----------------------------
//用Remez交错算法求解逼近问题
void remez1(); // <-- grid[ngrid], des[ngrid], wt[ngrid], iext[nfcns+1]
remez1();
//----------------------------
//计算理想响应函数 [2] Fig.2(Continued)-----P.5
//[2] NEG==0, 300:
if (nodd==0){//goto 310
//310: 偶数
// [1]
h_[0] =0.25*alpha[nfcns]; //下标为 [0~nfilt-1]
h_[nfilt-1] =h_[0];
for (j=1; j<nm1; j++){
h_[j] =0.25*(alpha[nfcns-j]+alpha[nz-j]);
h_[nfilt-1-j] =h_[j];
}//315
h_[nm1] =0.5*alpha[1] +0.25*alpha[2]; //?待验证
h_[nfcns] =h_[nm1];
/*[2] 对比 (半边)
h[1] =0.25*alpha[nfcns];
for (j=2; j<=nm1; j++){
h[j] =0.25*(alpha[nfcns+2-j]+alpha[nz-j]);
}//315
h[nfcns] =0.5*alpha[1] +0.25*alpha[2];
*/
}else{
//[1]
for (j=0; j<nm1; j++){//下标为 [0~nfilt-1]
h_[j] =0.5*alpha[nfcns-j];
h_[nfilt-1-j] =h_[j];
}//305
h_[nm1] =alpha[1];
