关于傅立叶算法的一个细节问题(快速,蝶形)。在线等待。

Kingore 2004-01-14 09:57:14
以下是傅立叶算法:(用在图象变换中)
在// 采用蝶形算法进行快速付立叶变换 这部分我看不懂:
向高手请教:
我把这部分代码COPY上来,如下:
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-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;
}

在蝶形算法中,我看过算法的大致原理:主要是:一个求N点的DFT可以被转换成两个
求N/2点的DFT。一直这样分解(都是蝶形形式);我们这里N=8;
最后一步是(因为是逆推,所以相当于算法中的第一步):
A[0]=X[0]+W[0]*X[4];
A[5]=X[0]-W[0]*X[4];

A[1]=X[1]+W[1]*X[5];
A[6]=X[1]-W[1]*X[5];
。。。。

而在上述算法程序中,到最内层的循环可以看出(设r=3,此时bfsize=4;p=0...)
所以可以得到最初的表达式为:
X2[0]=X1[0]+X1[4];
X2[4]=(X1[0]-X1[4])*W[0]

X2[1]=X1[1]+X1[5];
X2[5]=(X1[1]-X1[5])*W[1];
。。。
跟上面的分析并不一样!!!!
为什么?
各位大哥给我分析一下原因,谢谢。

(是我理解错了?)
还是使用了什么技巧?


以下是完整的程序:
///////////////////////////////////////////////////////////////////////
VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
// 付立叶变换点数
LONG count;

// 循环变量
int i,j,k;

// 中间变量
int bfsize,p;

// 角度
double angle;

complex<double> *W,*X1,*X2,*X;

// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
W = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-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];
}

// 释放内存
delete W;
delete X1;
delete X2;
}
BOOL WINAPI Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{

// 指向源图像的指针
unsigned char* lpSrc;

// 中间变量
double dTemp;

// 循环变量
LONG i;
LONG j;

// 进行付立叶变换的宽度和高度(2的整数次方)
LONG w;
LONG h;

int wp;
int hp;

// 图像每行的字节数
LONG lLineBytes;

// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行付立叶变换的宽度和高度(2的整数次方)
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
}

// 分配内存
complex<double> *TD = new complex<double>[w * h];
complex<double> *FD = new complex<double>[w * h];

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 指向DIB第i行,第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 给时域赋值
TD[j + w * i] = complex<double>(*(lpSrc), 0);
}
}

for(i = 0; i < h; i++)
{
// 对y方向进行快速付立叶变换
FFT(&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方向进行快速付立叶变换
FFT(&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;
}

// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行,第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j,是为了将变换后的原点移到中心
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}

// 删除临时变量
delete TD;
delete FD;

// 返回
return TRUE;
}

...全文
116 点赞 收藏 8
写回复
8 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复
Kingore 2004-01-14
我已将帖子加分,现在是100分!!。
各位兄弟帮忙,我郁闷死了。
为什么?

以下是傅立叶算法:(用在图象变换中)
在// 采用蝶形算法进行快速付立叶变换 这部分我看不懂:
向高手请教:
我把这部分代码COPY上来,如下:
// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-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;
}

在蝶形算法中,我看过算法的大致原理:主要是:一个求N点的DFT可以被转换成两个
求N/2点的DFT。一直这样分解(都是蝶形形式);我们这里N=8;
最后一步是(因为是逆推,所以相当于算法中的第一步):
A[0]=X[0]+W[0]*X[4];
A[5]=X[0]-W[0]*X[4];

A[1]=X[1]+W[1]*X[5];
A[6]=X[1]-W[1]*X[5];
。。。。

而在上述算法程序中,到最内层的循环可以看出(设r=3,此时bfsize=4;p=0...)
所以可以得到最初的表达式为:
X2[0]=X1[0]+X1[4];
X2[4]=(X1[0]-X1[4])*W[0]

X2[1]=X1[1]+X1[5];
X2[5]=(X1[1]-X1[5])*W[1];
。。。
跟上面的分析并不一样!!!!
如果象上面分析那样的话,应该为:
X2[0]=X1[0]+X1[4]*W[0];
X2[4]=X1[0]-X1[4]*W[0];

为什么?
各位大哥给我分析一下原因,谢谢。

(是我理解错了?)
还是使用了什么技巧?
回复
panadaa 2004-01-14
顶顶顶,关注。
回复
Kingore 2004-01-14
我要疯拉!,顶!!
回复
Kingore 2004-01-14
没有人耐心看看吗?
做过这方面的应该看看就明白了吧?
回复
Kingore 2004-01-14
大家不要看到代码长就放弃,只要看前面一部分就可以了,
后面的代码是因为完整性我才帖上的:)多谢!
回复
yaolan1999 2004-01-14
顶,这么长。
回复
Kingore 2004-01-14
自己顶
回复
panadaa 2004-01-14
⒉ 频域分组(按频率抽取的FFT算法)
回复
发帖
图形处理/算法
创建于2007-09-28

1.9w+

社区成员

VC/MFC 图形处理/算法
申请成为版主
帖子事件
创建了帖子
2004-01-14 09:57
社区公告
暂无公告