在优快云关于FIR的内容已经很多了，但值得收藏的寥寥无几。
本人感觉代码在此似乎大受欢迎，文献资料并未受到重视。所以我建议注重应用效率的朋友，比如只爱Matlab的朋友，可不必关注以下内容。
前一段时间研究信号分析书本上介绍的FIR设计最佳逼近法，内容都点到了，看得似懂非懂，待编写程序时发现很多细节不可回避，书中难以面面俱到。虽然也提供了程序，但缺少说明的Fortran代码就如同孔明布下的八卦阵，入者不死也要大伤元气。
前几天有幸在这里(http://liuxiaoyi125762092.download.youkuaiyun.com/)下载到最好的资料:
H.McCLELLAN, A computer program for designing optimum FIR linear phase digital filters.
经过几天发狠调试，总算把程序调通了，但仍没有时间来全部消化这个Fortran程序的全部细节和奥妙。现在放到这里，有兴趣的不妨来研究研究，交流心得。
/*
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]中也可能有错误, 如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];
/*[2] 对比 (半边) OK
for (j=1; j<=nm1; j++){
h[j] =0.5*alpha[nz-j];
}//305
h[nfcns] =alpha[1];
*/
}
/* NEG==1
if (nodd==0){//goto L330;
h[1] =0.25*alpha[nfcns];
for (j=2; j<=nm1; j++){
h[j] =0.25*(alpha[nz-j]-alpha[nfcns+2-j]);
}//335
h[nfcns] =0.5*alpha[1] -0.25*alpha[2];
}else{//(nodd!=0)
h[1] =0.25*alpha[nfcns];
h[2] =0.25*alpha[nm1];
for (j=3; j<nm1; j++){
int nzmj =nfcns -j;
h[j] =0.25*(alpha[nz-j]-alpha[nfcns+3-j]);
}//325
h[nfcns] =0.5*alpha[1] -0.25*alpha[3];
h[nz] =0;
//goto L350;
}
*/
vdaa_putFloat(h+1, nfilt, chView-1, 0, "h", "");
printf("/n/t Deviation Deviation(db)/n");
for (j=1; j<=nbands; j++){
float a =dev/wtx[j];
printf("/nBand%d %f %f", j, a, 20.*log10(a+fx[j]));
}//370
printf("/n/n/t滤波器系数:/n");
for (j=1; j<=nfcns; j++)
printf("/nh[%d] %f", j-1, h[j]);
//nfilt =nfilt +1; [1]
return ;
}
float gee(float y[], float x[], double ad[], int n, float grid[], int k);
float geeA(float y[], float x[], double ad[], int n, float f);
//
void remez1()
//
{
//以下注释主要来自文献[2]的程序流程图
float a[66+1], p[66+1], q[66+1];
int itrmax =25; //假设最大迭代次数
int niter =0;
float devl =-1.; //偏差
int nz =nfcns +1;
int nzz =nfcns +2;
float err; //局部临时变量
int j, jchnge;
double y1=0;
double ynz=0, comp=0;
int luck, kup, k, k1, nut1, nut, nu, l;
int knz, klow;
double dnum =0.;
double dden =0.;
double aa, bb, dtemp;
L100:
iext[nzz] =ngrid +1;
if (itrmax<=niter++) goto L400;
//计算Lagrange插值点横坐标 x[j] =cos(2Pi*F(j)), F(j)是极值点频率
for (j=1; j<=nz; j++){
x[j] =cos(grid[iext[j]]*pi2); //grid[iext[j]] 是极值点频率
}//110
//计算Lagrange插值系数(用函数d), ad[j] = 连乘(x[i]-x[j]), i=1~n且i!=j
{
int jet =(nfcns -1)/15 +1; //jet为适当的步长以抽取若干插值点来估计偏差? 这可能是一个技巧, 避免连乘的数据量过大导致上溢或下溢
double d(float x[], int k, int n, int m);
for (j=1; j<=nz; j++){
ad[j] =d(x, j, nz, jet);
}//120
}
//计算当前值的偏差 dev (即 [1] P.253 (8.4.8a) lo)
dnum =0.;
dden =0.;
k =1;
for (j=1; j<=nz; j++){
int i =iext[j];
dnum +=ad[j]*des[i];
dden +=k*ad[j]/wt[i];
k =-k;
} //130
dev =dnum/dden;
//Record sgn(dev) & set dev =|dev|
nu =(0.<dev)? -1: 1; //循环初值
if (dev<0.) dev =-dev; //dev *=-nu;
//计算Ragrange插值点纵坐标 y[j] = DES[iext[j]]+ (-1)^j dev/wt[iext[j]]
k =nu;
for (j=1; j<=nz; j++){
k =-k;
int i =iext[j];
//y[j] =des[i] +k*dev/wt[i]; [2],[1]源码, 应属笔误
//2项之间加号应为减号 -, 见 [3] P.528 (8.277), P.527 (8.2.72); [1] P.253 (8.4.8d). (否则测试失败)
y[j] =des[i] -k*dev/wt[i];
//实际上系数k在第1次循环为 sgn(dev), 然后每次反转——道理又何在? 待测试相反的起始符号
}//140
if (debugView){
put2Vdaa(iext+1, iext+1, nfcns+1, ngrid, chView+1, "iext", ""); //本次迭代初始的极值点，也是上次迭代的结果
//Pause! 显示之前的误差曲线"Error"和之前的极值点位置(曲线"iext_0"),
// 观察本次迭代如何确定了新的极值点(曲线"iext")
put2Vdaa(iext+1, iext+1, nfcns+1, ngrid, chView+0, "iext_0", ""); //上次的极值点位置
vdaa_putFloat(x+1, nz, chView+2, 0, "x", "");
vdaa_put1D(ad+1, nz, chView+3, 0, "ad", "");
put2Vdaa(y+1, iext+1, nz, ngrid, chView+4, "y", "");
}
{//先准备好所有格点上的偏差值
for (int i=1; i<=ngrid; i++){
float yp =gee(y, x, ad, nz, grid, i);
Error[i] =(yp-des[i])*wt[i];
Hr[i] =yp;
}
for (j=1; j<=nz; j++){//节点
int i =iext[j];
Hr[i] =y[j];
Error[i] =(y[j]-des[i])*wt[i];
}
if (debugView){
vdaa_putFloat(Error+1, ngrid, chView+5, 0, "Error", ""); //需搜索极值点的误差曲线
vdaa_putFloat(Hr+1, ngrid, chView+6, 0, "Hr", "");
//以上观察点[1~ngrid]未包括过渡带, 故重新获得全貌如下
float g =0, dg =0.5/(ngrid-1);
float *HrA =Hr;
int i=0;
while (i<ngrid){
HrA[i] =geeA(y, x, ad, nz, g); //此时HrA[i]与Hr[i]下标无关
//但过渡带的偏差不确定, 因未给定理想响应中过渡带的值
g +=dg;
i++;
}
vdaa_putFloat(HrA, ngrid, chView+7, 0, "HrA", "");
}
}
//dev < dev last time ? yes-> OUCH error message -> 计算最佳逼近系数;
//else
// No -> For the 1st extremal frequency set upper KUP and lower KLOW limits on the Grid points to be searched KUP =IEXT[2] KLOW = 0
//Furthermore, let the sign of the error at each extremal frequency be denoted by cigma(j) --?
//if (devl<dev) goto L150;
if (dev<devl){
//Call ouch; void ouch(){ printf("出错..."); }
printf("/n出错! 可能是累计摄入误差引起.");
goto L400; //迭代结束
}
//L150:
//搜索n+1个极值点频率
devl =dev; //本次偏差(最大极值)供下次比较
jchnge =0; //变动了的局部极值点数量
k1 =iext[1];
knz =iext[nz];
klow =0;
nut =-nu;
j =1; //自第1个极值点开始
//SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION
//
//用gee()估计在第j个极值点的临近格点grid[iext[j]+1]的频率响应
//然后计算加权误差: err = wt[k+1]*(gee(k+1)-des[k+1]), 其中 k =iext[j]
//这是迭代中最麻烦处，也是本程序对Remez算法的发挥。正宗的Remez算法一次只改变1个交错点(最大的极值点)，迭代速度很慢。
L200:
//以下代码与P.7 Fig.5(continued)一致性较差, 应重新改写! 首先尽可能简化 P.7 Fig.5.
//尤其nut很不直观(大概是P.7 Fig.5中cigma[j]), nut是如何与极值的符号对应的, 待考察.
// 按后面的应用, nut应是偏差的符号, 既如此直接取符号则可读性更好.
//再者, 取局部极值功能可抽出为子函数. 但要解决1个问题——这个极值点对应于哪个新的交错点iext[j] ?
//comp的应用跨度太大，交叉太多
//关于求err的2行可用Error[L]取代
if (j==nzz) ynz =comp;
if (nzz<=j) goto L300; //对j结束循环
kup =iext[j+1];
l =iext[j] +1;
nut =-nut;
if (j==2) y1 =comp;
comp =dev;
/* [2]
if (kup<=l) goto L220;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (nut*err<=comp) goto L220; // |err| < |dev| ?
comp =nut*err; //?是否与绝对值等价 nut == sign(err), nut*err == |err|?
*/ if (nut*Error[l]<=dev || kup<=l) goto L220; //iext[j]右边附近没有所希望的极值
comp =nut*Error[l];
L210:
//iext[j]右边有大的差值
/* [2]
l++;
if (kup<=l) goto L215;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (nut*err<=comp) goto L215;
comp =nut*err;
goto L210;
*/
{ double &errlast =comp;
while (l++<kup){
double err =nut*Error[l];
if (err<errlast) break;
errlast =err;
}//找到局部极大值点l-1, 或没找到则 l=kup
//可以改写为顺绝对值递增前进的子函数，去掉这个上蹿下跳的comp
}
L215:
iext[j] =l-1; //更新极值点为[l-1]
j++;
klow =l-1;
jchnge++;
goto L200;
L220:
l--;
L225:
l--;
if (l<=klow) goto L250;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (comp<nut*err) goto L230;
if (jchnge<=0.) goto L225;
goto L260;
L230:
comp =nut*err;
L235:
/* [2]
l--;
if (l<=klow) goto L240;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (nut*err-comp<=0.) goto L240;
comp =nut*err;
goto L235;
*/
{ double &errlast =comp;
while (klow<l--){
double err =nut*Error[l];
if (err<errlast) break;
errlast =err;
}
}
L240:
klow =iext[j];
iext[j] =l+1; //更新极值点为[l+1]
j++;
jchnge++;
goto L200;
L250:
l =iext[j] +1;
if (0.<jchnge) goto L215;
L255:
l++;
if (kup<=l) goto L260;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (nut*err-comp<=0.) goto L255;
comp =nut*err;
goto L210;
L260:
klow =iext[j];
j++;
goto L200;
L300:
if (nzz<j) goto L320;
if (iext[1]<k1) k1 =iext[1];
if (knz<iext[nz]) knz =iext[nz];
nut1 =nut;
nut =-nu;
l =0;
kup =k1;
comp =ynz*1.00001; //大概为了避免对走过的路径判断不一致?
luck =1;
L310:
l++;
if (kup<=l) goto L315;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (nut*err<=comp) goto L310;
comp =nut*err;
j =nzz;
goto L210;
L315:
luck =6;
goto L325;
L320:
if (9<luck) goto L350;
if (y1<comp) y1 =comp;
k1 =iext[nzz];
L325:
l =ngrid +1;
klow =knz;
nut =-nut1;
comp =y1*1.00001;
L330:
l--;
if (l<=klow) goto L340;
err =gee(y, x, ad, nz, grid, l);
err =(err-des[l])*wt[l];
if (nut*err<=comp) goto L330;
j =nzz;
comp =nut*err;
luck +=10;
goto L235;
L340:
if (luck==6) goto L370;
for (j=1; j<=nfcns; j++){
iext[nzz-j] =iext[nz-j];
}//345
iext[1] =k1; //?
goto L100;
L350:
{
int kn =iext[nzz];
for (j=1; j<=nfcns; j++){
iext[j] =iext[j+1];
}//360
iext[nz] =kn;
}
goto L100;
L370:
if (0.<jchnge) goto L100;
//CALCULATION OF THE COEFFIClENTS OF THE BEST APPROXIMATION
//USING THE INVERSE DISCRETE FOURIER TRANSFORM
// 用DFT计算最佳逼近系数alpha[]
L400:
int nm1 =nfcns -1;
float fsh =1.e-06;
float gtemp =grid[1];
x[nzz] =-2.;
float cn =2*nfcns -1;
float delf =1./cn;
l =1;
int kkk =0;
//if (edge[1]==0.0 && 0.5==edge[2*nbands]) kkk =1; -- [2]
if (grid[1]<=0.01 && 0.49<grid[ngrid]) kkk =1;
if (nfcns <= 3) kkk =1;
if (kkk == 1) goto L405;
dtemp =cos(pi2*grid[1]);
dnum =cos(pi2*grid[ngrid]);
aa =2./(dtemp -dnum);
bb =-(dtemp +dnum)/(dtemp -dnum);
L405:
for (j=1; j<=nfcns; j++){//->430
float ft =(j-1)*delf;
float xt =cos(pi2*ft);
if (kkk == 1) goto L410;
xt =(xt-bb)/aa;
ft =atan2(sqrt(1.-xt*xt), xt)/pi2;
L410:
float xe =x[l];
if (xe<xt) goto L420;
if ((xe-xt)<fsh) goto L415;
l++;
goto L410;
L415:
a[j] =y[l];
goto L425;
L420:
if ((xt-xe)<fsh) goto L415;
grid[1] =ft;
a[j] =gee(y, x, ad, nz, grid, 1);
L425:
if (1<l) l--;
L430: continue;
}
grid[1] =gtemp;
dden =pi2/cn;
for (j=1; j<=nfcns; j++){ //->510
dtemp =0;
dnum =(j-1)*dden;
if (nm1<1) goto L505;
for (k=1; k<=nm1; k++){
dtemp +=a[k+1]*cos(k*dnum);
}//500
L505: dtemp =2*dtemp +a[1];
L510: alpha[j] =dtemp;
}//510
for (j=2; j<=nfcns; j++)
alpha[j] *=2./cn; //550
alpha[1] /=cn;
if (kkk == 1) goto L545;
p[1] =2.*alpha[nfcns]*bb +alpha[nm1];
p[2] =2.*alpha[nfcns]*aa;
q[1] =alpha[nfcns-2] -alpha[nfcns];
for (j=2; j<=nm1; j++){//540
if (j<nm1) goto L515;
aa *=.5;
bb *=.5;
L515:
p[j+1] =0;
for (k=1; k<=j; k++){
a[k] =p[k];
p[k] =2.*bb*a[k];
}//520
p[2] +=a[1]*2.*aa;
for (k=1; k<j; k++){
p[k] +=q[k] +aa*a[k+1];
}//525
for (k=3; k<=j+1; k++){
p[k] +=aa*a[k-1];
}//530
if (j == nm1) goto L540;
for (k=1; k<=j; k++){
q[k] =-a[k];
}//535
q[1] +=alpha[nfcns-1-j];
L540: continue;
}//540
for (j=1; j<=nfcns; j++){
alpha[j] =p[j];
}//543
L545:
if (3<nfcns) return ;
alpha[nfcns+1] =0.0;
alpha[nfcns+2] =0.0;
return ;
}
//计算Lagrange插值系数(供gee()使用)
// 本函数与 [2] P.6(/21) Fig.5 (或 [1] P.253 (8.4.8b)——应去除(-1)^k) 并不相同!
double d(float x[], int k, int n, int m)
{
double d_ =1.;
float q =x[k];
for (int L=1; L<=m; L++){
for (int j=L; j<=n; j+=m){ //Frotran 中的 Do 循环语句与 Basic 中的For类似: DO 循环末行标号[,] 循环控制变量=初值, 终值[, 增量]
if (j==k) continue;
d_ *=2.*(q-x[j]);
}
}
/*
//!假设实际应为
for (int j=1; j<=n; j++){
if (j==k) continue;
d_ *=(double)(q-x[j]); //(x[j]-q)与(q-x[j])当n为奇数时结果相差负号. (q-x[j])是正宗写法
}
*/
d_ =1./d_;
return d_;
}
//任意f处插值结果,f={0.~0.5}
float geeA(float y[], float x[], double ad[], int n, float f)
{
//插值的定义域实际上是 x ={cos(0~Pi), 其中(0~Pi)为不包括端点的开区间内格点集}
double xf =cos(pi2*f);
double p =0;
double d =0;
for (int j=1; j<=n; j++){
if (Equal(xf,x[j])) continue;
double c =ad[j]/(xf-x[j]);
d +=c;
p +=c*y[j];
}
return p/d;
}
//第k点(频率grid[k])处插值结果
float gee(float y[], float x[], double ad[], int n, float grid[], int k)
{
//return geeA(x, x, ad, n, grid[k]);
//插值的定义域实际上是 x ={cos(0~Pi), 其中(0~Pi)为不包括端点的开区间内格点集}
double xf =cos(pi2*grid[k]);
double p =0;
double d =0;
for (int j=1; j<=n; j++){
if (iext[j]==k) continue; //根据定义 x[j]=cos(grid[iext[j]]*pi2), 很明显要使 (xf-x[j])!=0 则必须 iext[j]!=k
//本主程序中不会出现这种条件(故[2]无此行), 但作为通用函数应避免这种可能(误对节点计算插值)
double c =ad[j]/(xf-x[j]);
d +=c;
p +=c*y[j];
}
return p/d;
}
//测试 [1] P.258, 260
int main()
{
vdaa_ready();
vdaa_setAutoRefresh();
const float PI =3.1415926;
//P.260 例 8.4.2
//对工频及其倍频陷波：50,100,150Hz
//设采样率 fs = 500Hz, 3个阻带的归一化频率中心为 0.1, 0.2, 0.3
int nbands =7;
float edge[] ={0.0,0.07, 0.09,0.11, 0.13,0.17, 0.19,0.21, 0.23,0.27, 0.29,0.31, 0.33,0.5};
float fx[] ={1.0, .0, 1.0, .0, 1.0, .0, 1.0};
float wtx[] ={ 8, 1, 8, 1, 8, 1, 8};
const int nfilt =65;
float h[nfilt];
defir3(nfilt, nbands, edge, fx, wtx, h);
return 0;
}
扫一扫
7828
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
