宁波建设局网站,做水果代理去哪个网站,外网门户网站建设方案,做电影网站还能赚钱吗设计思想 傅里叶变换#xff0c;表示能将满足一定条件的某个函数表示成三角函数#xff08;正弦和/或余弦函数#xff09;或者它们的积分的线性组合。在不同的研究领域#xff0c;傅里叶变换具有多种不同的变体形式#xff0c;如连续傅里叶变换和离散傅里叶变换。 快速傅…设计思想 傅里叶变换表示能将满足一定条件的某个函数表示成三角函数正弦和/或余弦函数或者它们的积分的线性组合。在不同的研究领域傅里叶变换具有多种不同的变体形式如连续傅里叶变换和离散傅里叶变换。 快速傅里叶变换 (Fast Fourier Transform)即利用计算机计算离散傅里叶变换DFT的高效、快速计算方法的统称简称FFT于1965年由J.W.库利和T.W.图基提出。 对多项式 f ( x ) ∑ i 0 n a i x i , g ( x ) ∑ i 0 n b i x i f(x)∑_{i0}^na_i x_i ,g(x)∑_{i0}^nb_i x_i f(x)∑i0n​ai​xi​,g(x)∑i0n​bi​xi​定义其乘积fg为 ( f g ) ( x ) ( ∑ i 0 n a i x i ) ( ∑ i 0 n b i x i ) (fg)(x)(∑_{i0}^na_i x_i )(∑_{i0}^nb_i x_i) (fg)(x)(∑i0n​ai​xi​)(∑i0n​bi​xi​) 显然我们可以以 O ( n 2 ) O(n^2) O(n2)的复杂度计算这个乘积的每一项的系数。 但FFT可以以 O ( n l o g n ) O(nlogn) O(nlogn)的时间复杂度来计算这个乘积。 快速傅立叶算法核心的思想是分治即把一个复杂的问题分解为一个小的类似问题进行求解。 假定待变换离散时间序列信号长度为 n 2 m n2^m n2m将X(n)按照奇偶分组 X ( k ) ∑ r 0 N / 2 − 1 x ( 2 r ) W N 2 r k ∑ r 0 N / 2 − 1 x ( 2 r 1 ) W N 2 r 1 k X(k)\sum_{r0}^{N/2-1}x(2r) W_N^2rk∑_{r0}^{N/2-1}x(2r1) W_N^{2r1}k X(k)∑r0N/2−1​x(2r)WN2​rk∑r0N/2−1​x(2r1)WN2r1​k 上式可变换为 X ( k ) ∑ r 0 N / 2 − 1 x ( 2 r ) W N / 2 r k W N k ∑ r 0 N / 2 − 1 x ( 2 r 1 ) W N / 2 r k X(k)∑_{r0}^{N/2-1}x(2r) W_{N/2}^{rk} W_N^k ∑_{r0}^{N/2-1}x(2r1) W_{N/2}^{rk} X(k)∑r0N/2−1​x(2r)WN/2rk​WNk​∑r0N/2−1​x(2r1)WN/2rk​ 令 A ( k ) ∑ r 0 N / 2 − 1 x ( 2 r ) W N / 2 r k A(k)∑_{r0}^{N/2-1}x(2r) W_{N/2}^{rk} A(k)∑r0N/2−1​x(2r)WN/2rk​ B ( k ) ∑ r 0 N / 2 − 1 x ( 2 r 1 ) W N / 2 r k B(k)∑_{r0}^{N/2-1}x(2r1) W_{N/2}^{rk} B(k)∑r0N/2−1​x(2r1)WN/2rk​ 其中 k 取 0 , 1 , … … N / 2 − 1 k取0,1,……N/2-1 k取0,1,……N/2−1,从而 X ( k ) A ( k ) W N k B ( k ) X(k)A(k)W_N^k B(k) X(k)A(k)WNk​B(k) 由于 A ( k ) , B ( k ) A(k),B(k) A(k),B(k)都是 N / 2 N/2 N/2点的DFTX(k)为N点的DFT这一分治思想还可以进一步做下去。这一方法通常是使用递归实现的并行优化的难度较高。 为并行化快速傅里叶变换需要使用非迭代版本即先预处理每个位置上元素变换后的位置每个位置分治后的最终位置为其二进制翻转后得到的位置然后先将所有元素移到变换后的位置之后直接循环合并。 调整完循环顺序之后第一层循环变量i表示每一层变换的跨度第二层循环变量j表示每一层变换的第一个起点第三层循环遍历k则表示实际变换的位置k和ki。在这里从第二层开始是没有循环依赖的即对于不同的j算法不会对同一块地址进行访问因为访问的下标kj(mod i)。 为公平起见用于对比的串行版本快速傅里叶变换是直接在并行版本上删去编译推导#pragma omp for得到的。这是因为递归版本的快速傅里叶变换通常有较大的函数递归开销。 主函数运行时传入三个参数第一个参数为.exe文件第二个参数为并行部分使用的线程数量第三个参数为傅里叶变换的幂次因为傅里叶变换算法本身要求长度为2的幂次。 调整线程数与幂次考虑并行化的加速比。 运行结果 幂次20 线程数23456串行时间(s)0.5350.5450.4750.5370.503并行时间(s)0.3250.2960.2650.2330.241加速比1.6461.8411.7922.3052.087 结果分析 观察结果发现对于并行程序虽然线程数增加但加速比始终保持在2-3之间猜想可能是以下原因OpenMP的parallel region结束时线程之间需要同步有隐式路障。最好的OpenMP使用场景是线程之间没有很多需要锁保护的共享访问parallel region应该尽可能大以抵消OpenMP多线程带来的额外同步开销。 程序源码 #include iostream #includevector #includecmath #includecomplex.h #include omp.h using namespace std; typedef long long ll; typedef double lf; #define M_PI 3.14159265358979323846 struct Rader : vectorint {Rader(int n) : vectorint(1 int(ceil(log2(n)))){for (int i at(0) 0; i size(); i)if (at(i) at(i 1) 1, i 1)at(i) size() 1;} }; struct FFT : Rader {vectorcomplexlf w;FFT(int n) : Rader(n), w(size(), polar(1.0, 2 * M_PI / size())){w[0] 1;for (int i 1; i size(); i)w[i] * w[i - 1];}vectorcomplexlf pfft1(const vectorcomplexlf a) const{vectorcomplexlf x(size()); #pragma omp parallel forfor (int i 0; i a.size(); i)x[at(i)] a[i];for (int i 1; i size(); i 1) #pragma omp parallel forfor (int j 0; j i; j)for (int k j; k size(); k i 1){complexlf t w[size() / (i 1) * j] * x[k i];x[k i] x[k] - t, x[k] t;}return x;}vectorcomplexlf pfft2(const vectorcomplexlf a) const{vectorcomplexlf x(size()); #pragma omp parallel #pragma omp forfor (int i 0; i a.size(); i)x[at(i)] a[i];for (int i 1; i size(); i 1) #pragma omp forfor (int j 0; j i; j)for (int k j; k size(); k i 1){complexlf t w[size() / (i 1) * j] * x[k i];x[k i] x[k] - t, x[k] t;}return x;}vectorcomplexlf fft(const vectorcomplexlf a) const{vectorcomplexlf x(size());for (int i 0; i a.size(); i)x[at(i)] a[i];for (int i 1; i size(); i 1)for (int j 0; j i; j)for (int k j; k size(); k i 1){complexlf t w[size() / (i 1) * j] * x[k i];x[k i] x[k] - t, x[k] t;}return x;}vectorll ask(const vectorll a, const vectorll b) const{vectorcomplexlf xa(a.begin(), a.end()), xb(b.begin(), b.end());xa fft(xa), xb fft(xb);for (int i 0; i size(); i)xa[i] * xb[i];vectorll ans(size());xa fft(xa), ans[0] xa[0].real() / size() 0.5;for (int i 1; i size(); i)ans[i] xa[size() - i].real() / size() 0.5;return ans;} }; int main(int argc, char **argv) {//if (argc 3)//return cerr Error: No Enough parameters ( argv[0] num-of-threads power-of-transform-length).\n, 0;//omp_set_num_threads(atoi(argv[1]));//FFT fft(1 atoi(argv[2]));FFT fft(1 20);for(int i2;i11;i){omp_set_num_threads(i);cout#iendl;vectorcomplexlf a(fft.begin(), fft.end());double t0 omp_get_wtime();vectorcomplexlf b fft.fft(a);double t1 omp_get_wtime();cout Serial Time: t1 - t0 s\n;vectorcomplexlf c fft.pfft1(a);double t2 omp_get_wtime();cout Parallel Time1: t2 - t1 s\n;vectorcomplexlf d fft.pfft2(a);double t3 omp_get_wtime();cout Parallel Time2: t3 - t2 s\n;if (b ! c||c!d)cerr Error: Parallel result are not equivalent to Serial result.\n;} } 参考资料 快速傅里叶变换FFT超详解 手把手教快速傅立叶变换FFT算法 OpenMP实现并行快速傅里叶变换