我想在你这个代码的基础上实现福利叶变换,用的是C++中的一段代码,但是不成功啊,很郁闷,你能看看怎么改吗
void ImageProcess::fft(QImage Image)
{
// 中间变量
double dTemp;
// 循环变量
int i;
int j;
// 进行付立叶变换的宽度和高度(2的整数次方)
int w;
int h;
int wp;
int hp;
//int byteWidth;
//BYTE *pImageDataTemp=(BYTE *)malloc(sizeof(BYTE)*biDataSize);
w = 1;
h= 1;
wp = 0;
hp = 0;
//std::complex<double> *TD = (std::complex<double>*)malloc(sizeof(std::complex<double>)*biDataSize);
//std::complex<double> *FD = (std::complex<double>*)malloc(sizeof(std::complex<double>)*biDataSize);
/*byteWidth = biWidth*3;
if (byteWidth%4)
byteWidth += 4-(byteWidth%4);*/
while(w * 2 <=biWidth)
{
w *= 2;
wp++;
}
while(h * 2 <= biHeight)
{
h *= 2;
hp++;
}
// biDataSize=w*h;
// 分配内存
std::complex<double> *TD = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(w*h));
std::complex<double> *FD = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(w*h));
BYTE *pByte=(BYTE *)malloc(sizeof(BYTE)*(w*h));
// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// if(8 == depth) //采用了256色调色板,8位颜色索引
//{
pByte = Image.bits() +(biHeight-1- i )* lLineBytes + j;
TD[j + w * i] = std::complex<double>(*(pByte), 0);
//}
/*else if(32 == depth)//32位表示,数据格式为0xFFBBGGRR或0xAABBGGRR
{
pByte = Image.bits() + i * lLineBytes + j * 4;
//根据RGB模式转化成YIQ色彩模式的方式,取Y作为灰度值
BYTE pixelValue = (BYTE)(0.299 * (float)pByte[0] + 0.587 * (float)pByte[1]
+ 0.114 * (float)pByte[2]);
TD[j + w * i] = std::complex<double>( pixelValue, 0);
}
else
{
return;
}*/
//memcpy(Image.bits(), pByte, biDataSize);
//}
}
for(i = 0; i < h; i++)
{
// 对y方向进行快速付立叶变换
onfft(&TD[w * i], &FD[w * i], wp);
}
// 保存变换结果
for(i = 0; i < h; i++)
{
for(j = 0; j <w; j++)
{
TD[i + h * j] = FD[j + w * i];
}
}
for(i = 0; i <w; i++)
{
// 对x方向进行快速付立叶变换
onfft(&TD[i * h], &FD[i * h], hp);
}
// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 计算频谱
dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() +
FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的,直接设置为255
dTemp = 255;
}
*pByte = dTemp;
/* if(8 == depth) //采用了256色调色板,8位颜色索引
{
pByte = Image.bits() + i * lLineBytes + j;
*pByte = dTemp;
}
else if(32 == depth)//32位表示,数据格式为0xFFBBGGRR或0xAABBGGRR
{
pByte = Image.bits() + i * lLineBytes + j * 4;
pByte[0] = pByte[1] = pByte[2] = dTemp;
}
else
{
return;
}*/
// 更新源图像
memcpy(Image.bits(), pByte, sizeof(BYTE)*(w*h));
}
}
// 删除临时变量
free (TD);
free (FD);
free(pByte);
}
}
void ImageProcess::onfft(std::complex<double> * TD, std::complex<double> * FD, int r)
{
/*if(0 == power)
{
y[0] = a[0];
return;
}
int n = 1 << power;
double angle = 2 * PI / n;
std::complex<double> wn(cos(angle), sin(angle));
std::complex<double> w(1, 0);
std::complex<double> *a0 = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(n/2));
std::complex<double> *a1 = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(n/2));
std::complex<double> *y0 = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(n/2));
std::complex<double> *y1 = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(n/2));
for(int i = 0; i < n / 2; i ++)
{
a0 = a[2 * i];
a1 = a[2 * i + 1];
}
//分开成两个子fft过程
onfft(a0, y0, power - 1);
onfft(a1, y1, power - 1);
std::complex<double> u;
for(int k = 0; k < n / 2; k++) //蝶形算法
{
u = w * y1[k];
y[k] = y0[k] + u;
y[k + n / 2] = y0[k] - u;
w = w * wn;
}
free( a0);
free(a1);
free( y0);
free( y1);*/
// 付立叶变换点数
LONG count;
// 循环变量
int i,j,k;
// 中间变量
int bfsize,p;
// 角度
double angle;
std::complex<double> *W,*X1,*X2,*X;
// 计算付立叶变换点数
count = 1 << r;
// 分配运算所需存储器
W =(std::complex<double>*)malloc(sizeof(std::complex<double>)*(count/2));
X1 = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(count));
X2 = (std::complex<double>*)malloc(sizeof(std::complex<double>)*(count));
// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W =std:: complex<double> (cos(angle), sin(angle));
}
// 将时域点写入X1
memcpy(X1, TD, sizeof(std::complex<double>) * count);
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)//第k次迭代
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);//第k次蝶型单元数目
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}
// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}
// 释放内存
free( W);
free (X1);
free( X2);