提高正弦/余弦和大型阵列的计算速度
Posted
技术标签:
【中文标题】提高正弦/余弦和大型阵列的计算速度【英文标题】：Improving computation speed for sine/cosine and large arrays 【发布时间】：2016-02-19 19:39:02 【问题描述】：对于信号处理，我需要计算相对较大的 C 数组，如下面的代码部分所示。到目前为止，这工作正常，不幸的是，实施速度很慢。 “calibdata”的大小约为 150k，需要针对不同的频率/相位进行计算。有没有办法显着提高速度？在 MATLAB 中对逻辑索引执行相同的操作要快得多。
我已经尝试过的：
使用正弦的泰勒近似：没有显着改善。 使用 std::vector，也没有显着的改进。代码：
double phase_func(double* calibdata, long size, double* freqscale, double fs, double phase, int currentcarrier)
for (int i = 0; i < size; i++)
result += calibdata[i] * cos((2 * PI*freqscale[currentcarrier] * i / fs) + (phase*(PI / 180) - (PI / 2)));
result = fabs(result / size);
return result;
最好的问候， 托马斯
【问题讨论】：
为了确定，您是否尝试在编译器选项中打开优化？ 也许可以离线计算并使用查找表（将 CPU 换成内存！） 本题请选择C和C++之一。 除了性能问题，最好初始化result。
您的准确度要求是什么？
【参考方案1】：
在优化代码以提高速度时，第 1 步是启用编译器优化。我希望你已经这样做了。
第 2 步是分析代码并准确查看时间是如何花费的。如果没有分析，您只是在猜测，最终可能会尝试优化错误的东西。
例如，您的猜测似乎是cos 函数是瓶颈。但另一种可能是角度的计算是瓶颈。以下是我将如何重构代码以减少计算角度所花费的时间。
double phase_func(double* calibdata, long size, double* freqscale, double fs, double phase, int currentcarrier)
double result = 0;
double angle = phase * (PI / 180) - (PI / 2);
double delta = 2 * PI * freqscale[currentcarrier] / fs;
for (int i = 0; i < size; i++)
result += calibdata[i] * cos( angle );
angle += delta;
return fabs(result / size);
【讨论】：
我在想类似的事情，我怀疑使用angle += delta而不是计算i*...来累积数值误差。
@Bob__ 这是一个很好的观点。鉴于size 为150K，可以粗略估计错误。假设每次加法都会在angle 中引入半个 LSB 的错误，那么在 150K 次加法之后，错误累积到大约 17 个 LSB。由于 double 具有 52 位精度，因此错误使 35 位正确。因此，误差是 340 亿分之一，这对于大多数信号处理任务来说已经足够了。【参考方案2】：
好吧，我可能会因为这个答案而受到鞭挞，但我会为此使用 GPU。因为您的数组似乎不是自引用的，所以对于大型数组，您将获得的最佳加速是通过并行化......到目前为止。我不使用 MATLAB，但我只是在 MathWorks 网站上快速搜索了 GPU 利用率：
http://www.mathworks.com/company/newsletters/articles/gpu-programming-in-matlab.html?requestedDomain=www.mathworks.com
在 MATLAB 之外，您可以自己使用 OpenCL 或 CUDA。
【讨论】：
老实说，我认为这是最好的答案。到目前为止，您可以使用编译器优化。我可能会建议 OP 首先尝试对代码进行多线程处理（减少可以很容易地进行多线程处理），并在转向 GPU 处理之前看看这是否会有所改善，特别是因为 GPU 处理通常保证 100-1000 的开销毫秒，无论所涉及的实际操作/数据大小如何。 @Xirema 毫秒是正确的吗？那是 0.1 到 1.0 秒，对吧？ @user3386109 根据我的经验，毫秒是正确的。当他们说“实施缓慢”时，OP 的含义并不明显。如果您正在处理网络/套接字数据，那么 100-1000 毫秒的延迟就不算什么了。如果您在嵌入式处理器上处理实时数据，则 100-1000 毫秒可以杀死您的应用程序。因此，对于 100-1000 毫秒代表显着延迟的不确定性，我正在对冲我的建议。 @Xirema 谢谢！在考虑将代码移动到 GPU 之前，这是一个很好的了解。 @Xirema 哎呀！我要提到为多线程分解数组，但是却陷入了寻找 MATLAB 的东西。接得好。另外，我知道有延迟，但这些数字太疯狂了！如果是这样的话，那么他们一定是在用 OpenGL 计算着色器做一些特殊的巫术；关闭 VSync，我可以在大型阵列 (+10K) 上获得 200+ fps (GTX 780)。我会尽快使用计算着色器运行一些测试。无论如何，设置要容易得多。【参考方案3】：你在执行时间的敌人是：
师 函数调用（包括循环中的隐式调用） 访问不同领域的数据 操作不同的指令你应该研究数据驱动编程和有效使用数据缓存。
部门
无论是硬件支持还是软件支持部门都需要很长时间。如果可能，则通过更改数字基数或从循环中排除（如果可能）来消除。
函数调用
最有效的执行方法是顺序。处理器为此进行了优化。分支可能需要处理器执行一些额外的计算（分支预测）或重新加载指令缓存/流水线。浪费时间（可能会花费在执行数据指令上）。
对此的优化是使用循环展开和小函数内联等技术。还可以通过简化表达式和使用布尔代数来减少分支的数量。
访问不同区域的数据 现代处理器经过优化，可以处理本地数据（一个区域中的数据）。一个例子是用数据加载内部缓存。具体来说，加载一个 cache line 与数据。例如，如果数组中的数据在一个位置，而余弦数据在另一个位置，这可能会导致重新加载数据缓存，再次浪费时间。
更好的解决方案是将所有数据连续放置或连续访问所有数据。与其对余弦表进行许多不连续的访问，不如按顺序查找一批余弦值（之间没有任何其他数据访问）。
不同的指令
现代处理器在处理一批类似指令方面效率更高。例如，模式 load, add, store 在执行所有加载，然后是所有添加，然后是所有存储时对块更有效。
总结
这是一个例子：
register double result = 0.0;
register unsigned int i = 0U;
for (i = 0; i < size; i += 2)
register double cos_angle1 = /* ... */;
register double cos_angle2 = /* ... */;
result += calibdata[i + 0] * cos_angle1;
result += calibdata[i + 1] * cos_angle2;
上面的循环展开，类似的操作是分组执行的。
尽管关键字register 可能已被弃用，但建议编译器使用专用寄存器（如果可能）。
【讨论】：
两个问题：编译器不应该能够自行展开这个循环吗？您是否有任何来源声称现代处理器使用一批类似的指令会更好。大约半年前完成了关于这个主题的大学课程，我对此有点怀疑。在您的示例中，批量执行所有加载/存储可能会受到内存带宽/缓存可用性的限制。只要您有足够的 ALU 单元可用，批量执行加法运算才会有效。 编译器可能会根据优化展开循环，但它不会对指令进行分组（除非通过特定的、专门的编译指示提供更多信息）。我在 ARM7 嵌入式系统上运行的波形平滑函数上执行了此操作。我使用示波器进行了测量，发现它将执行时间缩短了 100 微秒以上。 我还对 Cortex A8 进行了基准测试，其中数据被内联处理（数组是本地的）与另一个模块中的数据。内联数据具有令人难以置信的更快的性能时间。将数据组织得更符合数据缓存也节省了更多的执行时间。所有分析均使用示波器进行。 非常感谢。在 Wikipedia 上阅读它，似乎大多数 ARM7 系统都没有使用乱序执行，这可能是您的实际结果与我的理论假设之间存在差异的原因。 如果您认为答案有帮助，请点击复选标记。【参考方案4】：您可以尝试使用基于复指数的余弦定义：
j^2=-1.
存储exp((2 * PI*freqscale[currentcarrier] / fs)*j) 和exp(phase*j)。评估 cos(...) 然后在 for 循环中恢复到几个产品和添加，sin()、cos() 和 exp() 只被调用几次。
实现如下：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#define PI 3.141592653589
typedef struct cos_plan
double complex* expo;
int size;
cos_plan;
double phase_func(double* calibdata, long size, double* freqscale, double fs, double phase, int currentcarrier)
double result=0; //initialization
for (int i = 0; i < size; i++)
result += calibdata[i] * cos ( (2 * PI*freqscale[currentcarrier] * i / fs) + (phase*(PI / 180.) - (PI / 2.)) );
//printf("i %d cos %g\n",i,cos ( (2 * PI*freqscale[currentcarrier] * i / fs) + (phase*(PI / 180.) - (PI / 2.)) ));
result = fabs(result / size);
return result;
double phase_func2(double* calibdata, long size, double* freqscale, double fs, double phase, int currentcarrier, cos_plan* plan)
//first, let's compute the exponentials:
//double complex phaseexp=cos(phase*(PI / 180.) - (PI / 2.))+sin(phase*(PI / 180.) - (PI / 2.))*I;
//double complex phaseexpm=conj(phaseexp);
double phasesin=sin(phase*(PI / 180.) - (PI / 2.));
double phasecos=cos(phase*(PI / 180.) - (PI / 2.));
if (plan->size<size)
double complex *tmp=realloc(plan->expo,size*sizeof(double complex));
if(tmp==NULL)fprintf(stderr,"realloc failed\n");exit(1);
plan->expo=tmp;
plan->size=size;
plan->expo[0]=1;
//plan->expo[1]=exp(2 *I* PI*freqscale[currentcarrier]/fs);
plan->expo[1]=cos(2 * PI*freqscale[currentcarrier]/fs)+sin(2 * PI*freqscale[currentcarrier]/fs)*I;
//printf("%g %g\n",creall(plan->expo[1]),cimagl(plan->expo[1]));
for(int i=2;i<size;i++)
if(i%2==0)
plan->expo[i]=plan->expo[i/2]*plan->expo[i/2];
else
plan->expo[i]=plan->expo[i/2]*plan->expo[i/2+1];
//computing the result
double result=0; //initialization
for(int i=0;i<size;i++)
//double coss=0.5*creall(plan->expo[i]*phaseexp+conj(plan->expo[i])*phaseexpm);
double coss=creall(plan->expo[i])*phasecos-cimagl(plan->expo[i])*phasesin;
//printf("i %d cos %g\n",i,coss);
result+=calibdata[i] *coss;
result = fabs(result / size);
return result;
int main()
//the parameters
long n=100000000;
double* calibdata=malloc(n*sizeof(double));
if(calibdata==NULL)fprintf(stderr,"malloc failed\n");exit(1);
int freqnb=42;
double* freqscale=malloc(freqnb*sizeof(double));
if(freqscale==NULL)fprintf(stderr,"malloc failed\n");exit(1);
for (int i = 0; i < freqnb; i++)
freqscale[i]=i*i*0.007+i;
double fs=n;
double phase=0.05;
//populate calibdata
for (int i = 0; i < n; i++)
calibdata[i]=i/((double)n);
calibdata[i]=calibdata[i]*calibdata[i]-calibdata[i]+0.007/(calibdata[i]+3.0);
//call to sample code
clock_t t;
t = clock();
double res=phase_func(calibdata,n, freqscale, fs, phase, 13);
t = clock() - t;
printf("first call got %g in %g seconds.\n",res,((float)t)/CLOCKS_PER_SEC);
//initialize
cos_plan plan;
plan.expo=malloc(n*sizeof(double complex));
plan.size=n;
t = clock();
res=phase_func2(calibdata,n, freqscale, fs, phase, 13,&plan);
t = clock() - t;
printf("second call got %g in %g seconds.\n",res,((float)t)/CLOCKS_PER_SEC);
//cleaning
free(plan.expo);
free(calibdata);
free(freqscale);
return 0;
使用gcc main.c -o main -std=c99 -lm -Wall -O3 编译。使用您提供的代码，在我的计算机上使用 size=100000000 需要 8 秒，而 建议的解决方案的执行时间需要 1.5 秒...这并不令人印象深刻，但也不容忽视。
提出的解决方案不涉及在 for 循环中对 cos 或 sin 的任何调用。事实上，只有乘法和加法。瓶颈要么是内存带宽，要么是通过平方求幂中的测试和对内存的访问（很可能是第一个问题，因为我添加了一个额外的复数数组）。
c中的复数见：
How to work with complex numbers in C? Computing e^(-j) in C如果问题是内存带宽，则需要并行性……直接计算cos会更容易。如果freqscale[currentcarrier] / fs 是整数，则可以进行额外的简化。您的问题非常接近Discrete Cosine Transform 的计算，目前的技巧接近离散傅里叶变换，而 FFTW 库非常擅长计算这些变换。
请注意，由于失去意义，当前代码可能会产生不准确的结果：size 很大时，result 可能比 cos(...)*calibdata[] 大得多。使用部分总和可以解决问题。
【讨论】：
嗯...不错的理论，但具体如何实现呢？ 可能对你没有任何帮助。评估两个 exp() 表达式的成本可能大于仅评估 cos() 的成本 引入复数来加速代码？ @EugeneSh。 ：实现这一点需要使用pow() 来计算指数...由于所有i 都需要指数的幂，因此需要对pow() 进行有限数量的调用。问题是这可能会导致精度问题......我要试试......
复指数，不是吗？【参考方案5】：
消除- (PI / 2) 的简单触发标识。这也比尝试使用machine_PI 的减法更准确。当值接近 π/2 时，这一点很重要。
cosine(x - π/2) == -sine(x)
使用const 和restrict：好的编译器可以利用这些知识执行更多优化。 （另见@user3528438）
// double phase_func(double* calibdata, long size,
// double* freqscale, double fs, double phase, int currentcarrier)
double phase_func(const double* restrict calibdata, long size,
const double* restrict freqscale, double fs, double phase, int currentcarrier)
某些平台使用 float 与 double 执行更快的计算，但精度损失是可以容忍的。 YMMV。双向配置文件代码。
// result += calibdata[i] * cos(...
result += calibdata[i] * cosf(...
尽量减少重新计算。
double angle_delta = ...;
double angle_current = ...;
for (int i = 0; i < size; i++)
result += calibdata[i] * cos(angle_current);
angle_current += angle_delta;
不清楚为什么代码使用long size 和int currentcarrier。我希望使用相同的类型并使用类型size_t。这对于数组索引来说是惯用的。 @Daniel Jour
反向循环可以允许与 0 进行比较，而不是与变量进行比较。有时会获得适度的性能提升。
确保充分启用编译器优化。
大家一起
double phase_func2(const double* restrict calibdata, size_t size,
const double* restrict freqscale, double fs, double phase,
size_t currentcarrier)
double result = 0.0;
double angle_delta = 2.0 * PI * freqscale[currentcarrier] / fs;
double angle_current = angle_delta * (size - 1) + phase * (PI / 180);
size_t i = size;
while (i)
result -= calibdata[--i] * sinf(angle_current);
angle_current -= angle_delta;
result = fabs(result / size);
return result;
【讨论】：
函数值和带缩放周期的正弦/余弦值的乘积之和看起来像卷积。我敢打赌phase_func2 函数会针对phase 的每个值重复调用。作者可能正在解决一个 XY 问题，其中 Y 是phase_func2 优化问题，X 是一种小波变换。
@mbaitoff 怀疑你是对的。像往常一样，来自 OP 的信息越多，优化的可能性就越大。【参考方案6】：
利用您拥有的内核，无需借助 GPU，而是使用 OpenMP。使用 VS2015 进行测试，优化器将不变量从循环中提取出来。启用 AVX2 和 OpenMP。
double phase_func3(double* calibdata, const int size, const double* freqscale,
const double fs, const double phase, const size_t currentcarrier)
double result;
constexpr double PI = 3.141592653589;
#pragma omp parallel
#pragma omp for reduction(+: result)
for (int i = 0; i < size; ++i)
result += calibdata[i] *
cos( (2 * PI*freqscale[currentcarrier] * i / fs) + (phase*(PI / 180.0) - (PI / 2.0)));
result = fabs(result / size);
return result;
启用 AVX 的原始版本耗时：~1.4 秒 并添加 OpenMP 将其降低到：~0.51 秒。
两个 pragma 和一个编译器开关的回报相当不错。
【讨论】：
以上是关于提高正弦/余弦和大型阵列的计算速度的主要内容，如果未能解决你的问题，请参考以下文章