急。。哪位有巴特沃斯带通滤波的算法？谢谢～

[guopo]

一定送分相谢！请发到[e-mail]或贴出来。

万分感谢！！

[guopo]

up!

[changlele]

又发送了一次，再次查收一下吧。应该没问题吧。

[changlele]

既然楼主说没找到代码，我把代码贴出来

/*************************************************************************
*
* 函数名称：
* DIBWinnerFilter()
*
* 参数:
* CDib *pDib - 指向CDib类的指针
*
* 返回值:
* BOOL - 维纳滤波复原操作成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用来对DIB图像进行维纳滤波复原操作。
*
************************************************************************/
void ButterWorthHighPass(LPBYTE lpImage, int nWidth, int nHeight, int nRadius)
{
// 循环控制变量
int y ;
int x ;

double dTmpOne ;
double dTmpTwo ;

// ButterWorth 滤波系数
double H ;

// 傅立叶变换的宽度和高度（2的整数次幂）
int nTransWidth ;
int nTransHeight;

double dReal ;
double dImag ;
// 图象象素值
unsigned char unchValue;

// 指向时域数据的指针
complex<double> * pCTData ;
// 指向频域数据的指针
complex<double> * pCFData ;

// 计算进行傅立叶变换的点数 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;

// 计算进行傅立叶变换的点数 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;

// 分配内存
pCTData=new complex<double>[nTransWidth * nTransHeight];
pCFData=new complex<double>[nTransWidth * nTransHeight];

// 初始化
// 图像数据的宽和高不一定是2的整数次幂，所以pCTData
// 有一部分数据需要补0
for(y=0; y<nTransHeight; y++)
{
for(x=0; x<nTransWidth; x++)
{
pCTData[y*nTransWidth + x]=complex<double>(0,0);
}
}

// 把图像数据传给pCTData
for(y=0; y<nHeight; y++)
{
for(x=0; x<nWidth; x++)
{
unchValue = lpImage[y*nWidth +x];
pCTData[y*nTransWidth + x]=complex<double>(unchValue,0);
}
}

// 傅立叶正变换
DIBFFT_2D(pCTData, nWidth, nHeight, pCFData) ;

// 下面开始实施ButterWorth高通滤波
for(y=0; y<nTransHeight; y++)
{
for(x=0; x<nTransWidth; x++)
{

H = (double)(y*y+x*x) ;
H = (nRadius * nRadius) / H ;
H = 1/(1+H) ;
pCFData[y*nTransWidth + x]=complex<double>(H*(pCFData[y*nTransWidth + x].real()),
H*(pCFData[y*nTransWidth + x].imag()) );
}
}

// 经过ButterWorth高通滤波的图象进行反变换
IFFT_2D(pCFData, pCTData, nWidth, nHeight);

// 反变换的数据传给lpImage
for(y=0; y<nHeight; y++)
{
for(x=0; x<nWidth; x++)
{
dReal = pCTData[y*nTransWidth + x].real() ;
dImag = pCTData[y*nTransWidth + x].imag() ;
unchValue = (unsigned char)max(0,min(255,sqrt(dReal*dReal+dImag*dImag)+100 ));
lpImage[y*nWidth + x] = unchValue ;
}
}


// 释放内存
delete pCTData;
delete pCFData;
pCTData = NULL;
pCFData = NULL;

}

傅立叶变换代码
/*************************************************************************
*
* \函数名称：
* IFFT_2D()
*
* \输入参数:
* complex<double> * pCFData - 频域数据
* complex<double> * pCTData - 时域数据
* int nWidth - 图像数据宽度
* int nHeight - 图像数据高度
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶反变换。
*
************************************************************************
*/
VOID WINAPI IFFT_2D(complex<double> * pCFData, complex<double> * pCTData, int nWidth, int nHeight)
{
// 循环控制变量
int x;
int y;

// 临时变量
double dTmpOne;
double dTmpTwo;

// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;

// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;

// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;

// 分配工作需要的内存空间
complex<double> *pCWork= new complex<double>[nTransWidth * nTransHeight];

//临时变量
complex<double> *pCTmp ;

// 为了利用傅立叶正变换,可以把傅立叶频域的数据取共轭
// 然后直接利用正变换，输出结果就是傅立叶反变换结果的共轭
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTmp = &pCFData[nTransWidth * y + x] ;
pCWork[nTransWidth * y + x] = complex<double>( pCTmp->real() , -pCTmp->imag() );
}
}

// 调用傅立叶正变换
::DIBFFT_2D(pCWork, nWidth, nHeight, pCTData) ;

// 求时域点的共轭，求得最终结果
// 根据傅立叶变换原理，利用这样的方法求得的结果和实际的时域数据
// 相差一个系数
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTmp = &pCTData[nTransWidth * y + x] ;
pCTData[nTransWidth * y + x] =
complex<double>( pCTmp->real()/(nTransWidth*nTransHeight),
-pCTmp->imag()/(nTransWidth*nTransHeight) );
}
}
delete pCWork ;
pCWork = NULL ;
}

[changlele]

/*************************************************************************
*
* \函数名称：
* FFT_2D()
*
* \输入参数:
* complex<double> * pCTData - 图像数据
* int nWidth - 数据宽度
* int nHeight - 数据高度
* complex<double> * pCFData - 傅立叶变换后的结果
*
* \返回值:
* 无
*
* \说明:
* 二维傅立叶变换。
*
************************************************************************
*/
VOID WINAPI DIBFFT_2D(complex<double> * pCTData, int nWidth, int nHeight, complex<double> * pCFData)
{

// 循环控制变量
int x;
int y;

// 临时变量
double dTmpOne;
double dTmpTwo;

// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;

// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;

// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;

// x，y（行列）方向上的迭代次数
int nXLev;
int nYLev;

// 计算x，y（行列）方向上的迭代次数
nXLev = (int) ( log(nTransWidth)/log(2) + 0.5 );
nYLev = (int) ( log(nTransHeight)/log(2) + 0.5 );

for(y = 0; y < nTransHeight; y++)
{
// x方向进行快速傅立叶变换
FFT_1D(&pCTData[nTransWidth * y], &pCFData[nTransWidth * y], nXLev);
}

// pCFData中目前存储了pCTData经过行变换的结果
// 为了直接利用FFT_1D，需要把pCFData的二维数据转置，再一次利用FFT_1D进行
// 傅立叶行变换（实际上相当于对列进行傅立叶变换）
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransHeight * x + y] = pCFData[nTransWidth * y + x];
}
}

for(x = 0; x < nTransWidth; x++)
{
// 对x方向进行快速傅立叶变换，实际上相当于对原来的图象数据进行列方向的
// 傅立叶变换
FFT_1D(&pCTData[x * nTransHeight], &pCFData[x * nTransHeight], nYLev);
}

// pCFData中目前存储了pCTData经过二维傅立叶变换的结果，但是为了方便列方向
// 的傅立叶变换，对其进行了转置，现在把结果转置回来
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransWidth * y + x] = pCFData[nTransHeight * x + y];
}
}

memcpy(pCTData, pCFData, sizeof(complex<double>) * nTransHeight * nTransWidth );
}