4000分求购一算法！并不很难，回者有分！

Question

有两个数组A[N],B[N](N为元素个数），
A[i]*B[i]=C[i]
然后d=C[0]+C[1]+。+C[N]
用快速傅立叶来实现！

Answer 1

在分布进行一维傅立叶变换时，多采用"蝴蝶图"的快速算法（详见信号处理方面资料），其核心算法如下：int N=(int)pow(2,M); file://N：序列长度（2的整数次幂）
ReverseOrder(A,N); file://对空间序列进行倒序
for(int i=1;i<=M;i++){
　int b=(int)pow(2,(i-1));
　for(int j=0;j<=(b-1);j++) {
　　float p=(float)(pow(2,(M-i))*j*2.0*PI/(float)N);
　　for(int k=j;k<=(N-1);){
　　　float tr=(float)(A[k+b].Re*cos(p)+A[k+b].Im*sin(p)); file://计算复数运算A*U
　　　float ti=(float)(A[k+b].Im*cos(p)-A[k+b].Re*sin(p));
　　　A[k+b].Re=A[k].Re-tr; file://复数运算A-tr
　　　A[k+b].Im=A[k].Im-ti;
　　　A[k].Re+=tr; file://复数运算A+tr
　　　A[k].Im+=ti;
　　　k+=b*2; 
　　}
　}
} 　　傅立叶逆变换的同傅立叶变换比较相似，只是在计算exp[j2πvy/N]时同正变换有符号的区别，以及在计算完成后逆变换需要将值除以N，因此不难写出一维傅立叶逆变换的实现代码。在进行二维傅立叶变换将图像由空域变换到频域之前，首先需要通过补0的手段将点数非2的整数次幂的非正方型网格采样构造为一个长宽均为2的整数次幂的正方型网格：int WM=(int)(log(W)/log(2)+1.0f); file://计算图像宽应为2的多少次幂
int HM=(int)(log(H)/log(2)+1.0f); file://计算图像高应为2的多少次幂
WM=HM=max(WM,HM); file://取二者大值
int WN=(int)pow(2,WM); file://构造网格宽度
int HN=(int)pow(2,HM); file://构造网格高度
for{int i=0;i;for(int j=0;j 　　if(i 　　　U[i*WN*3+j].Re=D[i*W*3+j]; file://D为图像序列
　　U[i*WN*3+j].Im=0.0f;
　　}else file://缺位补0
　　U[i*WN*3+j].Re=U[i*WN*3+j].Im=0.0f; 
　}
} 　　预处理完毕后，可对构造网格的每一列分别进行一维快速傅立叶变换，并将结果存放在原位置，结果乘以N，完成第一步的转换，求得F(x,v)：for(i=0;i 　for(int j=0;j 　　UH[j].Re=U[j*WN*3+i].Re; 
　　UH[j].Im=U[j*WN*3+i].Im; 
　}
　DFT_FFT(UH,HM); file://对UH进行快速离散傅立叶变换
　for(j=0;j 　　U[j*WN*3+i].Re=HN*UH[j].Re; file://N=HN
　　U[j*WN*3+i].Im=HN*UH[j].Im; 
　}
} 　　随即对构造网格的每一行进行傅立叶变换，得到最终的变换结果F(u,v)：for(i=0;i 　for(int k=0;k<3;k++){ file://对24位位图的R、G、B三分量均各自进行变换
　　for(int j=0;j 　　　UW[j].Re=U[i*WN*3+j*3+k].Re;
　　　UW[j].Im=U[i*WN*3+j*3+k].Im; 
　　}
　　DFT_FFT(UW,WM); file://对UW序列进行快速离散傅立叶变换
　　for(j=0;j 　　　U[i*WN*3+j*3+k].Re=UW[j].Re;
　　　U[i*WN*3+j*3+k].Im=UW[j].Im; 
　　}
　}
} 

Answer 2

我也有相同的问题。每次用FFT卷积下来都是个频域的数组，逆变换后是循环卷积的数组。
但程序只需要第一个循环卷积的结果。为节省运算量又不想让计算机运算其他的结果（数组太大了，用奔4 2.4G算起来都无法忍受）。感谢搂主的提问。
我也正在寻找办法，帮你顶。

Answer 3

Intel Performance library Suit

Answer 4

学习一下了！ipl库中有些函数不好用啊！iplGetPixel的到数值类型不对啊！
32F深度的图像，用iplGetPixel得到的数据非常恐怖！那个老兄使用过，或有类似的情况，解决了么？如何解决！
深深郁闷中

Answer 5

//chapter_01 基础数学概述
基础数学包含众多的科目，现在从“研究对象”和“发展方向”两个角度对其进行初步归纳
//01.01 按研究对象分类
1.算术/初等代数/高等代数/数论
2.欧式几何/非欧几何/ 解析几何/代数几何/微分几何/ 射影几何/分形几何/拓扑几何
2.实变函数论/复变函数论/泛函分析
3.微积分学/常微分方程/偏微分方程
5.概率和数理统计/突变理论/模糊数学/数理逻辑
6.运筹学/数学物理学/计算数学
//01.02 按发展方向分类
1.高等代数，拓展“数”的传统观念。新建代数的术语(如：群、环、域、模、格)并深入研究它们的特性
2.泛函分析，函数的传统定义是指两个数集之间的一个对应关系，现代数学将其理解为两个任意集合(不只是数集)之间的某个对应关系
3.拓扑几何，集合的概念也持续地进化着。人们曾经忽略了一维集合的内部结构；但出现二维集合后，人们发现无法把“四色猜想”归纳转化到其它更基本的科目中去；出现三维集合后，也无法把“莫比乌斯曲面”等价到其它科目的既有知识中。于是数学不得不增设堂口，另开门户，拜之为拓扑学
4.其它已经基本成熟了的科目，如：初等代数、欧式几何、微积分、概率论和数理统计本质上说，按研究对象的分类就足以准确描述其概貌了，但数学体系的各个分支是相互引证的，而当前某些的数学领域仍有一些问题尚未得到确证
所以说，某个方向上的重大进步可能可以完全改变整个数学的面貌，所以我们也需要基本了解发展方向
//01.03 注解
四色猜想，是指每幅地图都可以用四种颜色着色，使得有共同边界的国家都被着上不同的颜色
莫比乌斯曲面，将长纸带的一端在“垂直带长的平面内”翻转180度后，与另一端粘合，该纸带将不再具有两个侧面，新的曲面称为莫比乌斯曲面

Answer 6

//chapter_02 高等代数
//02.01 群、环、域的基本概念
当我们不再祈祷上帝，却跟从迦利略历行实验主义的时候，就该预料到世界将变得不可想象。因为我们选择了推理放弃了想象，我们选择做逻辑上正确的事而不是想象正确的事
术语定义的一致化和形式化，正是数学生命力的源泉。所以代数为不惜离开生动的形象，开始枯燥却符合逻辑的推导，这就是初等代数到高等代数的发展，本节介绍基本高等代数几个概念代数运算，集合S×S到S的一个映射就叫做集合S上的一个二元运算
代数系统，带有运算的集合叫做代数系统或代数结构
半群，设·是集合S上的一个二元运算，若满足对任意的a,b,c ∈S都有(a·b)·c = a·(b·c)。则称<S,·>是一个半群
群，设半群<G,·>具有单位元e，既，对任意的x∈G，总存在x'∈G，使x·x' = x'·x = e，则称半群<G,·>是一个群
交换群，若群<G,·>满足a·b=b·a， 则称<G,·>是交换群
既然数学的目标是“一致性和形式化的推导”，那么它并不能给我们呈现更多的新奇的东西，它做的只是尽量准确地描述基本概念
唯一特别的地方是群具有某个对称性，它是指：
对任意的x∈G，总存在x'∈G，使x·x' = x'·x = e
然后其它的这些生疏的术语只是在描述一些被用得烂熟的东西：符合结合律的运算(以及其所有涉及到的数)叫“半群”，满足交换律的运算(以及其所有涉及到的数)叫“交换群”环，具有两个二元运算“＋”与“·”的集合R所构成的代数结构(R,＋,·)满足以下条件，则称为之为环：
  ①<R,＋>是交换群
  ②<R,·>是半群
  ③对任意a,b,c ∈R且a·(b＋c)=a·b＋b·c，都有(a＋b)·c=a·c＋b·c
交换环，环(R,＋,·)的·运算满足a·b=b·a， 则称(R,＋,·)是交换环
除环，环(R,＋,·)中<R,·>也是群，不只是半群，则称环(R,＋,·)为除环
域，设环(R,＋,·)满足以下条件，则称之为域：
  ①<S,·>是群，记群<S,·>中的单位元为0，<S－{0},·>是交换群
如果只要一句话就可以准确理解群，那么那把拿群的概念组合三次后，引申出的的概念能有多复杂呢？环就是两个二元操作运算＋、·，其中＋满足交换律，·满足结合律，·对＋满足分配律的代数系统
定义不算很麻烦，但就是这样一个三句话的结构，却是数学理论的一大支柱，吸引了许多智商两百以上的数学家的穷其一生去追求，直到现在却仍未达完善。只要数学还有尚未解决的问题，高等代数就没得完善，高等代数也还在持续地发展
//02.02 (例子)李群研究
高等代数现在已经堪称内容丰富了。如其中相当活跃的科目李群(是一个连续的群，是以数学家Marius Sophus Lie名字命名的)研究，李先生在类比空间直线簇到和球面的对应关系后建立的从一个群到另一个群的变换关系，并利用这个变换来研究连续群的特性李先生发明的群变换具有正变换和逆变换两种形式：在源群内的某个数x，经正变换后得到了目标群内的数X，X再经逆变换将得到还原为原来的数x
只要存在这种可逆性，那我们在源群中的所有判断或命题都可以等价于目标群内的某个判断或命题，于是我们可以只研究目标群内的X的特性，就可以完全了解源群的特性
//02.03 数论
数论想证明的是“任何一个偶数n(n≥4)是两个素数之和”“任何整数都能唯一表示成几个素数的乘积”“x^n + y ^n = z^n在n>2时没有满足xyz≠0的整数解”这样一些命题
数论的途径是研究整除性等特性，它在计算机数据加密等方面具有广泛运用其实我看数学就是为了看数论变换，后来才发现按内涵而言的话，数论变换将会是“代数系统的变换”，比如群变换，或者域之类的什么变换，它是“群、环、域...”系列的中内容
我们也可以有专门的整数代数系统间的变换，但数论变换本身与数论是八竿子都打不着的关系
或者说数论变换具有广义和狭义之分，狭义定义中代数结构是以某数为模的整数环之间的变换，它可以应用于根据数字采样信号量，我们总可以把信号量（或许乘以某个常量后）取整近似，然后再分析描述被采样对象的特性
狭义数论变换除具有一般变换的性质外，还有一些依赖于数论可以证明的特性[也有人按照数学的发展历史把数论变换认为是狭义的，然后把其它一些变换如“拉普拉斯变换”“傅利叶变换”认为是积分变换中的内容，但它们在逻辑上是相等的，我看实在没必要区别∫和∑之间的关系]

Answer 7

//chapter_03 数论变换
//03.01 数论变换
李群中的变换具有正变换和逆变换两种形式：在源群内的某个数x，经正变换后得到了目标群内的数X，X再经逆变换将得到还原为原来的数x
由于这种唯一对应性我们可以得出一个结论：在源代数系统中的所有判断或命题都可以等价于目标代数系统内的某个判断或命题。于是我们也可以只研究目标代数系统内的X的特性，就可以完整地描述源代数系统的特性
这个“唯一对应性”就是数论变换的依据，“研究目标代数系统来描述源代数系统的特性”则是数论变换的工作方式数论变换的基本形式是如下一个积分式有限积分，其中积分区间为(a,b)：
F(α)=∫f(t)K(t, α)dt
K(t, α)可被称为核函数，不同的变换的核是不同的，一种变换对应于一个核函数
F(α)可被称为像函数
f(t)可被称为像原函数
(a,b)可被称为f(t)的定义域也就是说，偶对函数f(t)代入前述的积分式后，得到的新函数就是f(t)的变换的结果
该变换将“在代数系统t中的函数f(t)”变换为“在代数系统α中的函数F(α)”
最大的意外是被变换的对象不是数而是一个函数，这也不难理解，我们考虑代数系统时认为运算是其中的一个组成部分，于是如果只把代数系统作为集合处理显然不太恰当，而运算则也会体现出“集合的特性”，那数论变换把运算作为研究对象当然是合理的
从公式可见我在求解F(α)的一个值时，引用到了f(t)和K(t, α)在(a,b)区间上的所有取值，F(α)的值对于f(t)而言将体现出一种全局的统计性的规律，所以不必在试图想象像函数与像原间的关系，只要牢记前述转换公式即可
//03.02 (例子)拉氏变换
拉普拉斯变换是一种数论变换，其定义的核为K(t, α) = e^(-st)自然对数的-st次幂，f(t)的定义域为[0,+∞)
L(s)=∫f(t) * e^(-st)dt    其中t∈[0,+∞)拉氏变换主要用在系统的分析中，我们可以用拉氏变换构造出系统的传递函数G(s)
G(s) = L( out(t) ) / L(  in(t) )
可证 out(t) = L'( G(s)*L(in(t)) )
其中in(t)为系统输入，out(t)为系统输出，t为时间
考虑一个物理系统：在水平直轨道上有两个小车用弹簧k2连接，左侧小车用弹簧k1与地连接，左侧小车质量为m1，右侧小车质量为m2，现对m1施加向右的力f，求左小车的位移
m1*y1'' + k1*y1 + k2*(y1-y2) = f
m2*y2'' + k2*y2 = k2*y1
虽然是两个串联弹簧应该是相当简单的物理系统，在物理应用中我们很容易列出具有多个二次导数的方程，我是无法解出这个方程组
但如果你问我这个问题，我可以解出答案，因为我会拉普拉氏变换，上述方程组的像函数相对简单，所以我能够把它变换后求解再逆变换先代入微分算子，可得
m1*p^2*y1 + k1*y1 + k2*y1 -k2*y2 = f
(m2*p^2 + k2)*y2 = k2*y1
后式代入前式消去y2，可得
y1 = f * (m2*p^2 + k2) / [m1*m2*p^4 + (m1*k2 + k1*m2 + k2*m2)*p^2 + k1*k2]
其中in(t) = f，out(t) = y1，则有
G(s) = (m2*s^2 + k2)/[m1*m2*s^4 + (m1*k2 + k1*m2 + k2*m2)*s^2 + k1*k2]
既然out(f) = L'( G(s)*L(f) ) = L'( G(s)*f )
现在要做的就是要把G(s)代入上式然后就去查公式(调用函数库!)吧其实对我来说查拉氏变换公式表和查微分方程的公式表并没有什么不同，但拉普拉斯证明变换后的像函数L(s)有特殊的优势比如如下一个问题：
若偶在一个电压系统中加一个输出反馈，偶把输出电压反相后附加到输入电压上，要求求解新的输出电压与输入电压之间的关系
拉普拉斯证明这样的情况下 G_new(s) = G_old(s) / (1 + G_old(s))
类似的传递函数基本组合有多中，他们的出现充分证明了拉氏变换的实用性
//03.03 (例子)模M的整数环间的变换
本小节只讨论整数的情况（整数有许多的应用场合，如位图的像素点就是整型数据），先介绍数论中的两个概念
模M的整数环，例如计算机的整数是有字长限制的，考虑单字节整数会发现它只有0～255可取值，如果你要计算250+11的运算的话你只能得到5。如果当运算的结果不在0～255范围内时，我们强制执行取模操作（通过加或减256把值调整到0～255内），用取模的结果代替结果，这样的代数体系就是模M的整数环
本原单位根，当(α^N) % M = 1时（即，在模M的整数环中α^N = 1时），称α模M的N阶本原单位根
费马数，任何具有2^b + 1 (其中b=2^t，t是自然数)数论变换中若f(t)为若它是定义域为整数的循环函数，我们将其称为离散函数，显然任何离散函数f(t)，都可以被表述为具有周期个分量的向量
离散数论变换将具有一些特殊的性质和优点，我们将研究它，为了方便研究我们将再附加一个约束：相关函数的定义域和值区域都是模M的整数环
记f(t)函数的周期为N，显然N<M，则
f(t)可表述为一个N个整形分量的向量X
K(t, α)可表述为N阶的方形矩阵T={α(i,j)}
F(α)可表述为F=T*X
模M整型变换的主要参数包括三个M、N、α，于是如何确定这三个参数成了“模M的整数环间的变换”研究的焦点
从M、N的定义来看他们显然受f(t)的影响非常大，而在确定了M、N后，本原单位根的值也只有确定的几个了
当前确定M、N的最优化的思路是费马变换，它取M为某个费马数，N取2^t当变换的核函数T的为具有以下形式时：
α(i,j) = a^(i*j) a为M的N阶单位根
记F=T*f、H=T*h，利用不同参数的恒等式的特性，我们可以证明：
 F*H = T*(f¤h)
这就是变换的卷积性
//03.04 (附录)卷积的概念
(略)

Answer 8

//chapter_04 傅利叶变换
//04.01 优点
傅利叶变换也是一种具有卷积性的数论变换，变换的核函数为复函数(也可只取其实数部分)，这可能会令人疑惑：如果可以用整数变换实现一般的数论运算，那为何还要学习复数变换呢？我也不知道，反正流行是有目共睹的，我看见了有人问傅氏变换问题但从未见人问过费马变换
也有可能是傅氏变换更具有物理意义才使它如此流行，如果把像原函数认为是描述物理量的直观变化的话，那么傅氏变换的像函数将对应与该物理量的按三角函数的分解后在某频率上的振幅，也就是认为像原函数是有许多个正／余弦波叠加而成“以这些波的频率为定义域振幅为值域的函数”，例如
如果当像原函数描述位移时，变换的像函数描述该运动的振动的辆，也就是说像函数在某点上的取值将是这个频率上的振动的振幅
如果当像原函数描述电流时，像函数在某点上的取值将是该频率的交流电的电流大小
这样傅氏变换就成了物理学家的一个必备，因为在研究振动必须知道物体在各个频率上的振幅，预测LC 回路时就得先把电流都分解成交流电甚至于有些物理学家认为仅有这种死板的“按三角函数规律的分解”是不够的，给傅氏变换的单位根加了个矫正系数使单位根也沿定义域的变化而变化，称为小波变换。任何能用傅氏变换的地方均可用小波变换。
所以傅氏变换是一个岔路口，如果你只是想利用数论变换计算图象卷积的话，你应该去了解费马变换；如果你是被他的物理意义所打动的话，务必学习小波变换；没任何理由留在岔路口止步张望(偶是前者，浮点卷积也可以变形为整形的卷积嘛)
//04.02 离散傅利叶变换
让我们来了解非常流行的傅氏变换吧，我也就只介绍计算机形式的离散傅氏变换，离散就意味着像原函数的定义域是几个孤立的点，意味着偶可以把变换写成F=T*X 的矩阵形式
像原函数f(t)可表述为一个N个整形分量的向量X
核函数K(t, α)可表述为N阶的方形矩阵T={T(i,j) = α^(i*j)}　α=e^(-2*PI*j/N) i∈[0,N) j∈[0,N)
像函数F(α)可表述为F=T*X
(与“模M的整数环间的变换”不同的是其中的分量是实数而不是整数)
其中e为自然对数的底
PI为圆周率
j为-1的平方根
//04.03 快速傅利叶变换原理
当然F也是一个N维向量，记α=e^(-2*PI*j/N)，则F 的第k 个分量
F_k = ∑f_m * W^(m*k)  k ∈[0,N)
若N = N1 * N2，我们设
m=N1*m2+m1  m1,k1 ∈[0,N1)
k=N1*k1+k2  m2,k2 ∈[0,N2)
W1=e^(-2*PI*j/N1)
W2=e^(-2*PI*j/N2)
可得
F_N2*k1+k2 = ∑W_2 ^ (m2*k2) ∑(f_N1*m2+m1 * W^(m1*k2)) * W1 ^(m1*k1)   m1∈[0,N1-1) m2∈[0,N2-1)
易得原F_k 的复杂度为N*N
F_N2*k1+k2 的复杂度为N1*N2(N1 + N2 * 1)
如果N=2^t，那偶可以取N1=N2＝2^(t-1)，然后可以对N1,N2继续进行分解，到Nx=2为止
那么观察可的其复杂度为(N/2)log2(N)
我们就把原算法从N*N 优化到了(N/2)log2(N)这就是快速傅利叶变换的思路(从分解过程来看，这种优化方式也适合于“模M的整数环间的变换”)
//04.03 Cooley Tukey 算法
在1965 年Cooley 和Tukey 根据该思路优化，联合提出一个算法，他们绘制了一个信号流图表述，
记N=2^b，则共有b次跳跃，第t次的跳跃中，将第t-1列的转化为t列，该跳跃将处理2^(b-t+1)个元素
在由这2^(b-t+1)个分量组成的向量中，每间隔2^(b-t)个位置的两个分量组成一个独立的运算，一个第t次的跳跃由2^(b-t)个运算组成
第t次的跳跃的第r1个运算的定义是
第一步：
　f(col=t, row=r1) =         f(col=t-1, row=r2) + f(col=t-1, row=r)
　f(col=t, row=r2) = W_t^r *(f(col=t-1, row=r2) - f(col=t-1, row=r))
　其中
　 W_t = e^(-2*PI*j/N1)
　 r2 = r1 + 2^(b-t)
第二步：
　若2^(b-t+1) > 2（即，b>t），将第奇数个分量都提到第偶数个分量前[其实Cooley 和Tukey 提出有两套信号流图，区别在于W^(m1*k2)对应与时域抽选和和频域抽选，我这里引用的都是频域抽选]
//04.04 Winograd 算法
分析快速傅利叶变换原理中过程，我们发现优化的原因是因式分解，就是把大的变换运算分解为两个小的变换运算，再将小的运算的结果叠加模拟出大的变换运算
遵循着这条思路，但是因式分解的思路不一样就可以得到不同的快速傅利叶变换在04.02的分析中介绍的因式分解是Cooley Tukey算法的思路，现在我再介绍一个素因子FFT
F_k = ∑f_m * W^(m*k)  k ∈[0,N)
若N = N1 * N2，且N1,N2的最大公约数为1
则我们可以找到以下两个n和k
n=(N1*n2+N2*n1)%N  n1,k1 ∈[0,N1)
k=(N1*k2+N2*k1)%N  n2,k2 ∈[0,N2)
W1=e^(-2*PI*j/N1)
W2=e^(-2*PI*j/N2)
可得
F_N1*k2+N2*k1 = ∑∑f_N1*n2+N2*n1 * W1^(N2*n1*k1) * W2^(N1*n2*k2)   n1∈[0,N1-1) n2∈[0,N2-1)
然后，利用t2*k1及t1*k2来代替k1及k2，并使N2*t2%N1 = 1及N1*t1%N2 = 1，得
F_N1*t1*k2+N2*t2*k1 = ∑∑f_N1*n2+N2*n1 * W1^(n1*k1) * W2^(n2*k2)   n1∈[0,N1-1) n2∈[0,N2-1)该公式的意义是将一维变换转换为二维变换，据说该说法可以是Cooley Tukey 算法的五分之一

Answer 9

#define SWAP_NUMBER(type, A, B)  {type tTemp;tTemp=A;A=B;B=tTemp;}
#define SWAP_FLOAT(A, B)         SWAP_NUMBER(float,A,B)
static bool FFT_CooleyTukey (bool bForwordOrInvert, int nBounceMax, float *pfReal, float *pfImage)
{
  const long nCount = 1<<nBounceMax;
  {
    // Do the bit reversal 
    long j = 0;
    for (int i=0;i<nCount-1;i++)
    {
      if (i < j)
      {
        SWAP_FLOAT(pfReal [i], pfReal [j]);
        SWAP_FLOAT(pfImage[i], pfImage[j]);
      }      long k = nCount / 2;
      while (k <= j)
      {
        j -= k;
        k >>= 1;
      }
      j += k;
    }
  }  float  fSin  = -1.0;
  float  fCos  =  0.0;
  for (long nIndexBounce=0; nIndexBounce<nBounceMax; nIndexBounce++)
  {
    //bounce per loop
    {
      long  nIndex1 = 0==nIndexBounce ? 1 : 1<<nIndexBounce;
      long  nIndex2 = nIndex1 << 1;      float  u_1 = 1.0;
      float  u_2 = 0.0;
      for (long j=0;j<nIndex1;j++)
      {
        for (long nIndexTop=j; nIndexTop<nCount; nIndexTop+=nIndex2)
        {
          //an quadrangularly across operation
          const long nIndexDown = nIndexTop + nIndex1;
          const float        t1 = u_1 * pfReal [nIndexDown] - u_2 * pfImage[nIndexDown];
          const float        t2 = u_1 * pfImage[nIndexDown] + u_2 * pfReal [nIndexDown];
          pfReal [nIndexDown]   = pfReal [nIndexTop] - t1;
          pfImage[nIndexDown]   = pfImage[nIndexTop] - t2;
          pfReal [nIndexTop]    = pfReal [nIndexTop] + t1;
          pfImage[nIndexTop]    = pfImage[nIndexTop] + t2;
        }        float fTemp = u_1 * fSin - u_2 * fCos;
        u_2         = u_1 * fCos + u_2 * fSin;
        u_1         = fTemp;
      }
    }    fSin = float(sqrt((1.0 + fSin) / 2.0));
    fCos = float(sqrt((1.0 - fSin) / 2.0));
    if (bForwordOrInvert)
      fCos  *= -1;
  }  return true;
}

Answer 10

ANY PROBLEM？
偶可连续看了两周的书了