如何用CC++实现合成地震记录的完整流程？
- 内容介绍
- 文章标签
- 相关推荐
本文共计13469个文字，预计阅读时间需要54分钟。
C++实现地震记录合成+实例：从波阻抗模型中获取与地震对应的反射系数，再将反射系数与子波积分得到合成地震记录。+1+由波阻抗获取反射系数+地震波在介质中传播时，反射系数传递，作用于面。
C/C++实现合成地震记录本实例将从波阻抗模型中获得与之对应的反射系数，再将反射系数与子波褶积得到合成地震记录。
1 由波阻抗获取反射系数地震波在介质中传播时，作用于某个面积上的压力与单位时间内垂直通过此面积的质点流量(即面积乘质点振动速度）之比，具有阻力的含义，称为波阻抗，其数值等于介质密度\(\rho\)与波速\(v\)的乘积，即\(Z=\rho v\)。波阻抗差异形成地震反射，用反射系数表征上下界面差异，公式为
\[{r} = \frac{{{Z_2} - {Z_1}}}{{{Z_2} + {Z_1}}} = \frac{{{\rho _2}{v_2} - {\rho _1}{v_1}}}{{{\rho _2}{v_2} + {\rho _1}{v_1}}}\\Z_2=Z_1\frac{{1 + r}}{{1 - r}}\\ Z_{2}=Z_{1} \exp \left[\ln \left(\frac{1+r_{}}{1-r_{}}\right)\right] \]在实际情况下，反射系数\(r\)一般远远小于 1，即 \(r \ll 1\) 。则有如下关系：
\[\ln \left(\frac{1+r_{}}{1-r_{}}\right) \approx 2 r_{} \]上式可化为：
\[Z_{2} \approx Z_{1} e^{2 r} \\ r_1=\frac{1}{2}(\ln Z_2-\ln Z_1) \]于是，通过对上下两层波阻抗的对数值作差，就能由此得到对应的反射系数。
将上式用矩阵的形式表示为：
\[{\left[ {\begin{array}{*{20}{c}} { - \frac{1}{2}}&{\frac{1}{2}}&{}&{}&{}&{}\\ {}&{ - \frac{1}{2}}&{\frac{1}{2}}&{}&{}&{}\\ {}&{}& \ddots & \ddots &{}&{}\\ {}&{}&{}&{ - \frac{1}{2}}&{\frac{1}{2}}&{}\\ {}&{}&{}&{}&{ - \frac{1}{2}}&{\frac{1}{2}} \end{array}} \right]_{(n - 1) \times n}}\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\ln Z_1}}\\ {{\ln Z_2}} \end{array}}\\ \vdots \\ {{\ln Z_{n-1}}} \end{array}}\\ {{\ln Z_{n}}} \end{array}} \right]_{n\times 1} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{r_1}}\\ {{r_1}} \end{array}}\\ \vdots \\ {{r_{n-2}}} \end{array}}\\ {{r_{n-1}}} \end{array}} \right]_{(n-1)\times 1} \]容易看出整个实现过程，先是对波阻抗序列求对数，之后左乘差分矩阵，便能够得到反射系数矩阵（长度=\(n-1\)）
代码如下：
/*
输入：
波阻抗 float* Imp
波阻抗序列长度 int nt
输出：
反射系数 float* Rcs
*/
void Compute_Rcs(float *Imp, float* Rcs, int LenImp) {
float **Diff = NULL;
float *LnImp = NULL;
Diff = alloc2float(LenImp - 1, LenImp); //差分矩阵
LnImp = (float*)calloc(LenImp, sizeof(float)); //波阻抗对数值
for (int i = 0; i < LenImp - 1; i++) {
for (int j = 0; j < LenImp; j++) {
if (i == j) { //当矩阵行列下标相等(i=j)时，取值为-0.5
Diff[i][j] = -0.5;
}
else if (i == j - 1) { //当矩阵下标i=j-1时，取值为0.5
Diff[i][j] = 0.5;
}
else { //其他情况，值为0
Diff[i][j] = 0;
}
}
}
for (int i = 0; i < LenImp; i++) {
LnImp[i] = log(Imp[i]); //取对数
}
matrixMultilcation(Diff, LnImp, Rcs, LenImp, LenImp - 1); //矩阵乘法
}
2 褶积模型
地震道可以用褶积模型描述：
\[S(t) =w(t) *r(t) \]其中\(S(t)\)为地震道，\(w(t)\)为地震子波，\(r(t)\)为反射系数。震源在地表激发的时候震源在地表激发的时刻，地震子波被认为是一个尖脉冲信号。如果这个脉冲在传播过程中能量和波形均保持不变，那么地震记录就是反射系数剖面。而实际上子波在传播过程中，由于受球面扩散和吸收衰减的影响，并非严格的尖脉冲信号。
在本次实例中，使用的子波有以下三种：
/*
函数功能：
-由反射系数褶积地震子波，得到合成地震记录
输入：
-反射系数序列 float *Rcs
-反射系数长度 int LenRcs
-子波类型 int WaveletType
-子波长度 int WaveletLenTemp
-子波采样间隔 int WaveletInc
输出：
-合成地震记录 float *Syn
*/
void Rcs2Syn(float *Rcs, float *Syn, int LenRcs, int WaveletType, int WaveletLenTemp, int WaveletInc) {
float WaveletAmp = 1.0;
int WaveletLen = WaveletLenTemp /WaveletInc;
float *WaveletValue = NULL;
float *WaveletTime = NULL;
WaveletValue = (float*)calloc(WaveletLen, sizeof(float));
WaveletTime = (float*)calloc(WaveletLen, sizeof(float));
switch (WaveletType) {
case 1:
// type==1 雷克子波
parameterRicker Ricker_para;
Ricker_para.DomainFreq = 30;
Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f,
Ricker_para.DomainFreq, WaveletValue, WaveletTime, WaveletAmp);
break;
case 2:
//type==2 Yu's子波
parameterYus Yus_para;
Yus_para.QFreq = 60.0;
Yus_para.PFreq = 10.0;
Double_Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f, Yus_para.QFreq, Yus_para.PFreq,
WaveletValue, WaveletTime, WaveletAmp);
break;
case 3:
//type==3 双余弦子波
parameterDoubleCos DoubleCos_para;
DoubleCos_para.HighCut = 80;
DoubleCos_para.HighPass = 60;
DoubleCos_para.LowCut = 8;
DoubleCos_para.LowPass = 12;
int nK = getK(WaveletLen);
int SpecLen = (int)pow(2.0, nK - 1);
float FreqInc = 1000.0 / (2.0 * WaveletInc * SpecLen);
int ThickRatio = 8;
int DoubleCosSpecLen = int(SpecLen * 2 * ThickRatio);
float *DoubleCosFreq = NULL;
float *DoubleCosSpec = NULL;
DoubleCosFreq = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCosSpec = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCos_Spec(DoubleCosSpecLen, FreqInc / ThickRatio, DoubleCos_para.LowCut, DoubleCos_para.LowPass,
DoubleCos_para.HighPass, DoubleCos_para.HighCut, DoubleCosFreq, DoubleCosSpec, WaveletAmp);
DoubleCos_Wavelet(DoubleCosSpec, DoubleCosSpecLen, nK, ThickRatio, WaveletLen / 2, WaveletLen, WaveletInc / 1000.0,
WaveletTime, WaveletValue, 1.0);
free(DoubleCosFreq);
free(DoubleCosSpec);
break;
}
float *Signal = NULL;
int lenSignal = WaveletLen + LenRcs - 1;
Signal = (float*)calloc(lenSignal, sizeof(float));
convolve_cwp(LenRcs, 0, Rcs, WaveletLen, 0, WaveletValue, lenSignal, 0, Signal); //卷积
memcpy(Syn, Signal + (WaveletLen - 1) / 2, (LenRcs) * sizeof(float)); //赋值
free(WaveletValue);
free(WaveletTime);
free(Signal);
}
3 主函数
输入的波阻抗数据将被转换为合成记录：
int main(){
const char *filenameInput = "Imp.segy"; //输入segy文件
const char *filenameOutput = "Syn.segy"; //输出segy文件
bhed fileheader;
segy traceheader;
unsigned int nt = 0;
unsigned int sizefileheader = sizeof(fileheader);
unsigned int sizetraceheader = sizeof(traceheader);
unsigned int ncdp = 0;
long long size_file = 0;
long long size_trace = 0;
FILE *fpinput = NULL;
FILE *fpoutput = NULL;
float *dataInput = NULL;
float *dataOutput = NULL;
fpinput = fopen(filenameInput, "rb");
fpoutput = fopen(filenameOutput, "wb");
if (fpinput == NULL) {
printf("can not open %s file\n", filenameInput);
return false;
}
if (fpoutput == NULL) {
printf("can not open %s file\n", filenameOutput);
return false;
}
fread(&fileheader, sizefileheader, 1, fpinput);
nt = fileheader.hns;
_fseeki64(fpinput, 0, SEEK_END);
size_file = _ftelli64(fpinput);
size_trace = nt * sizeof(float) + sizetraceheader;
ncdp = (size_file - (long long)sizefileheader) / size_trace;
_fseeki64(fpinput, sizefileheader, SEEK_SET);
fwrite(&fileheader, sizefileheader, 1, fpoutput);
dataInput = (float*)calloc(nt, sizeof(float));
/****************以上为segy文件的读入*****************/
int WaveletLenTemp = 128; //子波长度
int WaveletInc = 1; //子波间隔
int WaveletType = 1; //选择子波类型，此处为Ricker
//遍历每一道
for (int itrace = 0; itrace < ncdp; itrace++)
{
fread(&traceheader, sizetraceheader, 1, fpinput);
fread(dataInput, nt * sizeof(float), 1, fpinput);
//预处理
int *start_end = NULL;
start_end = (int*)calloc(2, sizeof(int));
Preprocessing(dataInput, nt, start_end);//Preprocessing的作用是屏蔽0值和明显的异常值，读取有值的部分
int Start = 0;
int End = 0;
int LenImp = 0;
Start = start_end[0]; //获取有值部分的开始和结尾
End = start_end[1];
LenImp = End - Start + 1;
int Boundary = 3; //为了防止边界出现假轴，将波阻抗稍微延长出一部分
int LenImp_ex = LenImp + Boundary * 2;
float *Imp = NULL;
float *Rcs = NULL;
float *syn = NULL;
Imp = (float*)calloc(LenImp_ex , sizeof(float));
Rcs = (float*)calloc(LenImp_ex - 1, sizeof(float));
syn = (float*)calloc(LenImp_ex - 1, sizeof(float));
for (int i = Boundary; i < LenImp+Boundary; i++) {
Imp[i] = dataInput[Start-Boundary + i];
}
//上边界
for (int i = 0; i < Boundary; i++) {
Imp[i] = dataInput[Start];
}
//下边界
for (int i = LenImp+Boundary; i < LenImp + Boundary * 2; i++) {
Imp[i] = dataInput[End];
}
//第一步，计算反射系数
Compute_Rcs(Imp,Rcs,LenImp);
//第二步，反射系数褶积子波，得到合成记录
Rcs2Syn(Rcs,syn, LenImp_ex - 1, WaveletType, WaveletLenTemp, WaveletInc);
dataOutput = (float*)calloc(nt, sizeof(float));
for (int i = 0; i < LenImp; i++) {
dataOutput[Start + i] = syn[i + Boundary]; //传入合成记录
}
//写入segy
fwrite(&traceheader, sizetraceheader, 1, fpoutput);
fwrite(dataOutput, nt * sizeof(float), 1, fpoutput);
free(start_end);
free(Imp);
free(Rcs);
free(syn);
}
free(dataInput);
free(dataOutput);
fclose(fpinput);
fclose(fpoutput);
return 0;
}
完整代码
I fundamental.h
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include"complex.h"
#include"alloc.h"
void matrixMultilcation(float **data1, float *data2, float *dataOut, int c, int m); //矩阵乘法
void convolve_cwp(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z); //褶积
int getK(int N); //求2的K次方
void kkfft(float *pr, float *pi, int n, int k, float *fr, float *fi, int l, int il);//傅氏变换
II fundamental.cpp
#include "fundamental.h"
// 矩阵乘法
// data1 m行c列
//data2 c行1列
void matrixMultilcation(float **data1, float *data2, float *dataOut, int c, int m) {
float sumTemp = 0;
for (int i = 0; i < m; i++) {
sumTemp = 0;
for (int k = 0; k < c; k++) {
sumTemp += data1[i][k] * data2[k];
}
dataOut[i] = sumTemp;
}
}
// 满足2^k>=N的最小k值,返回k值
int getK(int N)
{
int k, temp = 1;
for (k = 1; k < N; k++)
{
temp *= 2;
if (temp >= N)
{
break;
}
}
return k;
}
void kkfft(float *pr, float *pi, int n, int k, float *fr, float *fi, int l, int il)
{
int it, m, is, i, j, nv, l0;
float p, q, s, vr, vi, poddr, poddi;
for (it = 0; it <= n - 1; it++)
{
m = it;
is = 0;
for (i = 0; i <= k - 1; i++)
{
j = m / 2;
is = 2 * is + (m - 2 * j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306 / (1.0 * n);
pr[1] = cos(p);
pi[1] = -sin(p);
if (l != 0)
{
pi[1] = -pi[1];
}
for (i = 2; i <= n - 1; i++)
{
p = pr[i - 1] * pr[1];
q = pi[i - 1] * pi[1];
s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
pr[i] = p - q;
pi[i] = s - p - q;
}
for (it = 0; it <= n - 2; it = it + 2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr + fr[it + 1];
fi[it] = vi + fi[it + 1];
fr[it + 1] = vr - fr[it + 1];
fi[it + 1] = vi - fi[it + 1];
}
m = n / 2;
nv = 2;
for (l0 = k - 2; l0 >= 0; l0--)
{
m = m / 2;
nv = 2 * nv;
for (it = 0; it <= (m - 1)*nv; it = it + nv)
{
for (j = 0; j <= (nv / 2) - 1; j++)
{
p = pr[m * j] * fr[it + j + nv / 2];
q = pi[m * j] * fi[it + j + nv / 2];
s = pr[m * j] + pi[m * j];
s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
poddr = p - q;
poddi = s - p - q;
fr[it + j + nv / 2] = fr[it + j] - poddr;
fi[it + j + nv / 2] = fi[it + j] - poddi;
fr[it + j] = fr[it + j] + poddr;
fi[it + j] = fi[it + j] + poddi;
}
}
}
if (l != 0)
{
for (i = 0; i <= n - 1; i++)
{
fr[i] = fr[i] / (1.0 * n);
fi[i] = fi[i] / (1.0 * n);
}
if (il != 0)
{
for (i = 0; i <= n - 1; i++)
{
pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]);
if (fabs(fr[i]) < 0.000001 * fabs(fi[i]))
{
if ((fi[i] * fr[i]) > 0)
{
pi[i] = 90.0;
}
else
{
pi[i] = -90.0;
}
}
else
{
pi[i] = atan(fi[i] / fr[i]) * 360.0 / 6.283185306;
}
}
}
}
return;
}
III convolution.cpp
/******************************************************************************
Input:
lx length of x array
ifx sample index of first x
x array[lx] to be convolved with y
ly length of y array
ify sample index of first y
y array[ly] with which x is to be convolved
lz length of z array
ifz sample index of first z
Output:
z array[lz] containing x convolved with y
******************************************************************************/
#include"fundamental.h"
/* prototype of function used internally */
static void convolve_cwp_s(int, int, float*, int, int, float*, int, int, float*);
void convolve_cwp(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z)
#ifdef SIMPLE_CONV
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1, i, j, jlow, jhigh;
float sum;
x -= ifx; y -= ify; z -= ifz;
for (i = ifz; i <= ilz; ++i) {
jlow = i - ily; if (jlow < ifx) jlow = ifx;
jhigh = i - ify; if (jhigh > ilx) jhigh = ilx;
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
}
#else
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1,
i, j, ilow, ihigh, jlow, jhigh;
float sa, sb, xa, xb, ya, yb, *t;
/* if x is longer than y, swap x and y */
if (lx > ly) {
i = ifx; ifx = ify; ify = i;
i = ilx; ilx = ily; ily = i;
i = lx; lx = ly; ly = i;
t = x; x = y; y = t;
}
/* handle short x with special code */
if (lx >= 1 && lx <= 30) {
convolve_cwp_s(lx, ifx, x, ly, ify, y, lz, ifz, z);
return;
}
/* adjust pointers for indices of first samples */
x -= ifx;
y -= ify;
z -= ifz;
/* OFF LEFT: i < ify+ifx */
/* zero output for all i */
ilow = ifz;
ihigh = ify + ifx - 1; if (ihigh > ilz) ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
/* ROLLING ON: ify+ifx <= i < ify+ilx */
/* if necessary, do one i so that number of j in overlap is odd */
if (i < ify + ilx && i <= ilz) {
jlow = ifx;
jhigh = i - ify;
if ((jhigh - jlow) % 2) {
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
}
/* loop over pairs of i and j */
ilow = i;
ihigh = ilx + ify - 1; if (ihigh > ilz) ihigh = ilz;
jlow = ifx;
jhigh = ilow - ify;
for (i = ilow; i < ihigh; i += 2, jhigh += 2) {
sa = sb = 0.0;
xb = x[jhigh + 1];
yb = 0.0;
for (j = jhigh; j >= jlow; j -= 2) {
sa += xb * yb;
ya = y[i - j];
sb += xb * ya;
xa = x[j];
sa += xa * ya;
yb = y[i + 1 - j];
sb += xa * yb;
xb = x[j - 1];
}
z[i] = sa;
z[i + 1] = sb;
}
/* if number of i is odd */
if (i == ihigh) {
jlow = ifx;
jhigh = i - ify;
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
/* MIDDLE: ify+ilx <= i <= ily+ifx */
/* determine limits for i and j */
ilow = i;
ihigh = ily + ifx; if (ihigh > ilz) ihigh = ilz;
jlow = ifx;
jhigh = ilx;
/* if number of j is even, do j in pairs with no leftover */
if ((jhigh - jlow) % 2) {
for (i = ilow; i < ihigh; i += 2) {
sa = sb = 0.0;
yb = y[i + 1 - jlow];
xa = x[jlow];
for (j = jlow; j < jhigh; j += 2) {
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
xa = x[j + 2];
}
z[i] = sa;
z[i + 1] = sb;
}
/* else, number of j is odd, so do j in pairs with leftover */
}
else {
for (i = ilow; i < ihigh; i += 2) {
sa = sb = 0.0;
yb = y[i + 1 - jlow];
xa = x[jlow];
for (j = jlow; j < jhigh; j += 2) {
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
xa = x[j + 2];
}
z[i] = sa + x[jhigh] * y[i - jhigh];
z[i + 1] = sb + x[jhigh] * y[i + 1 - jhigh];
}
}
/* if number of i is odd */
if (i == ihigh) {
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
/* ROLLING OFF: ily+ifx < i <= ily+ilx */
/* if necessary, do one i so that number of j in overlap is even */
if (i <= ily + ilx && i <= ilz) {
jlow = i - ily;
jhigh = ilx;
if (!((jhigh - jlow) % 2)) {
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
}
/* number of j is now even, so loop over both i and j in pairs */
ilow = i;
ihigh = ily + ilx; if (ihigh > ilz) ihigh = ilz;
jlow = ilow - ily;
jhigh = ilx - 2; /* Dave's new patch */
for (i = ilow; i < ihigh; i += 2, jlow += 2) {
sa = sb = 0.0;
xa = x[jlow];
yb = 0.0;
for (j = jlow; j < jhigh; j += 2) {
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
xa = x[j + 2];
}
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
z[i] = sa;
z[i + 1] = sb;
}
/* if number of i is odd */
if (i == ihigh) {
jlow = i - ily;
jhigh = ilx;
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
/* OFF RIGHT: ily+ilx < i */
/* zero output for all i */
ilow = i;
ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
}
/* internal function optimized for short x */
static void convolve_cwp_s(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z)
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1,
i, j, ilow, ihigh, jlow, jhigh;
float x0, x1, x2, x3, x4, x5, x6, x7, x8, x9,
x10, x11, x12, x13, x14, x15, x16, x17, x18, x19,
x20, x21, x22, x23, x24, x25, x26, x27, x28, x29,
ya, yb, z0, z1, sum;
x -= ifx;
y -= ify;
z -= ifz;
/* OFF LEFT: i < ifx+ify */
ilow = ifz;
ihigh = ify + ifx - 1; if (ihigh > ilz) ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
/* ROLLING ON: ify+ifx <= i < ify+ilx */
ilow = ify + ifx; if (ilow < ifz) ilow = ifz;
ihigh = ify + ilx - 1; if (ihigh > ilz) ihigh = ilz;
jlow = ifx;
jhigh = ilow - ify;
for (i = ilow; i <= ihigh; ++i, ++jhigh) {
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
/* MIDDLE: ify+ilx <= i <= ily+ifx */
ilow = ify + ilx; if (ilow < ifz) ilow = ifz;
ihigh = ily + ifx; if (ihigh > ilz) ihigh = ilz;
if (lx == 1) {
x0 = x[ifx];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 2) {
x0 = x[ifx];
x1 = x[ifx + 1];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 3) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 4) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 5) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 6) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 7) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 8) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 9) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 10) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 11) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 12) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 13) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 14) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 15) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 16) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 17) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 18) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 19) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 20) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 21) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 22) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 23) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 24) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 25) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 26) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 27) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 28) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
x27 = x[ifx + 27];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb; z1 += x27 * yb;
ya = y[i - ifx - 27]; z0 += x27 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 29) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
x27 = x[ifx + 27];
x28 = x[ifx + 28];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb; z1 += x27 * yb;
ya = y[i - ifx - 27]; z0 += x27 * ya; z1 += x28 * ya;
yb = y[i - ifx - 28]; z0 += x28 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 30) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
x27 = x[ifx + 27];
x28 = x[ifx + 28];
x29 = x[ifx + 29];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb; z1 += x27 * yb;
ya = y[i - ifx - 27]; z0 += x27 * ya; z1 += x28 * ya;
yb = y[i - ifx - 28]; z0 += x28 * yb; z1 += x29 * yb;
ya = y[i - ifx - 29]; z0 += x29 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
if (ihigh >= ilow && (ihigh - ilow) % 2 == 0) {
ilow = ihigh;
jlow = ifx;
jhigh = ilx;
for (i = ilow; i <= ihigh; ++i) {
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
}
/* ROLLING OFF: ily+ifx < i <= ily+ilx */
ilow = ily + ifx + 1; if (ilow < ifz) ilow = ifz;
ihigh = ily + ilx; if (ihigh > ilz) ihigh = ilz;
jlow = ilow - ily;
jhigh = ilx;
for (i = ilow; i <= ihigh; ++i, ++jlow) {
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
/* OFF RIGHT: ily+ilx < i */
ilow = ily + ilx + 1; if (ilow < ifz) ilow = ifz;
ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
}
#endif
#ifdef TEST
static void convolve_cwp1(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z)
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1, i, j, jlow, jhigh;
float sum;
x -= ifx; y -= ify; z -= ifz;
for (i = ifz; i <= ilz; ++i) {
jlow = i - ily; if (jlow < ifx) jlow = ifx;
jhigh = i - ify; if (jhigh > ilx) jhigh = ilx;
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
}
#include "cwp.h"
main()
{
int lx, ly, lz, ifx, ify, ifz, i;
float *x, *y, *z, *z1;
/* loop over tests */
while (1) {
/* make parameters */
lx = 1 + 100 * franuni();
ly = 1 + 100 * franuni();
lz = 1 + 100 * franuni();
ifx = -100 + 200 * franuni();
ify = -100 + 200 * franuni();
ifz = -100 + 200 * franuni();
/* allocate space */
x = alloc1float(lx);
y = alloc1float(ly);
z = alloc1float(lz);
z1 = alloc1float(lz);
/* fill x and y arrays */
for (i = 0; i < lx; i++)
x[i] = i + 1;
for (i = 0; i < ly; i++)
y[i] = i + 1;
/* convolve and check */
convolve_cwp1(lx, ifx, x, ly, ify, y, lz, ifz, z1);
convolve_cwp(lx, ifx, x, ly, ify, y, lz, ifz, z);
for (i = 0; i < lz; ++i)
if (z[i] != z1[i]) break;
if (i < lz) {
printf("%10d %10d %10d %10d %10d %10d\n",
lx, ifx, ly, ify, lz, ifz);
printf("z1[%d]=%g != z[%d]=%g\n",
i + ifz, z1[i], i + ifz, z[i]);
pp1d(stdout, "simple", lz, ifz, z1);
pp1d(stdout, "optimized", lz, ifz, z);
exit(-1);
}
/* free space */
free1float(x);
free1float(y);
free1float(z);
free1float(z1);
}
}
#endif
IV ImpToSyn.h
#include"fundamental.h"
using namespace std;
//雷克子波
void Ricker(int Lr, //中心点
int L, //子波长度
float dalt, //采样间隔，单位：s
float fp, //主频
float *w, //子波数组
float *time,//横轴数组
float amp); //最大振幅
void Double_Ricker(int Lr, int L, float dalt, float fq, float fp, float *w, float *time, float amp);//Yu'S
void DoubleCos_Spec(int L, float dalf, float fLowCut, float fLowPass, float fHighPass, float fHighCut, float *pFreq, float *pSpec, float amp);//双余弦
void DoubleCos_Wavelet(float *pSpec, int nFreqLen, int nK,int nThickRatio,int Lr,int L,float dalt,float *pTime, float *pValue, float amp);
void Compute_Rcs(float *Imp, float* Rcs,int LenImp); //计算反射系数
void Rcs2Syn(float *Rcs, float *Syn, int LenRcs, int WaveletType, int WaveletLenTemp,int WaveletInc); //反射系数褶积子波得到合成记录
void Preprocessing(float *Imp_in, unsigned int LenImp, int *start_end); //对输入地震道进行预处理
typedef struct parameterwavelet {
int type; // 子波类型
int WaveletLen;
int WaveletInc;
}parameterwavelet;
typedef struct parameterRicker { // 雷克子波参数
int DomainFreq;
}parameterRicker;
typedef struct parameterYus { // Yu's子波参数
float QFreq;
float PFreq;
}parameterYus;
typedef struct parameterDoubleCos { // DoubleCos子波参数
float LowCut;
float LowPass;
float HighPass;
float HighCut;
}parameterDoubleCos;
V ImpToSyn.cpp
#include"ImpToSynProfile.h"
#define Min_Imp 1000.0
/*雷克子波*/
void Ricker(int Lr, //中心点
int L, //子波长度
float dalt, //采样间隔，单位：s
float fp, //主频
float *w, //子波数组
float *time,//横轴数组
float amp) //最大振幅
{
double temp;
double pi;
int Ll;
int k, i;
Ll = L - Lr - 1;
pi = 4.0 * atan(1.0);
for (i = -Lr; i <= Ll; i++)
{
k = i + Lr;
temp = pi * fp * (dalt * i);
temp = temp * temp;
time[k] = (float)i * dalt * 1000.0;
w[k] = amp * (1.0 - 2.0 * temp) * exp(-temp);
}
return;
}
/*Yu's子波*/
void Double_Ricker(int Lr, int L, float dalt, float fq, float fp, float *w, float *time, float amp) //与Ricker相同，有Q,P两个频率
{
double temp1, temp2;
double pi;
int Ll;
int k, i;
Ll = L - Lr - 1;
pi = 4.0 * atan(1.0);
for (i = -Lr; i <= Ll; i++)
{
k = i + Lr;
temp1 = pi * fq * (dalt * i);
temp2 = pi * fp * (dalt * i);
temp1 = temp1 * temp1;
temp2 = temp2 * temp2;
time[k] = (float)i * dalt * 1000.0;
w[k] = amp * (1.0 / (fq - fp)) * (fq * exp(-temp1) - fp * exp(-temp2));
}
return;
}
//计算双余弦子波频谱
void DoubleCos_Spec(int L, float dalf, float fLowCut, float fLowPass, float fHighPass, float fHighCut, float *pFreq, float *pSpec, float amp)
{
//------ bak ------
int i, j;
float fqe_temp;
float PI;
float xishu = 0;
if (L % 2 > 0)
{
L--;
}
PI = 4.0 * atan(1.0);
for (i = 0; i < L / 2; i++)
{
fqe_temp = i * dalf;
pFreq[i] = fqe_temp;
if (fqe_temp <= fLowCut)
{
pSpec[i] = 0.0;
}
else if (fLowCut < fLowPass && fqe_temp > fLowCut && fqe_temp < fLowPass)
{
pSpec[i] = 1.0 * ((xishu + 1.0) / (1.0 - xishu)
+ cos((fqe_temp - fLowCut) * PI / (fLowPass - fLowCut) + PI))
/ ((xishu + 1.0) / (1.0 - xishu) + 1.0);
pSpec[i] *= amp;
}
else if (fqe_temp >= fLowPass && fqe_temp <= fHighPass)
{
pSpec[i] = 1.0;
pSpec[i] *= amp;
}
else if (fHighPass < fHighCut && fqe_temp > fHighPass && fqe_temp < fHighCut)
{
pSpec[i] = 1.0 * ((xishu + 1.0) / (1.0 - xishu)
+ cos((fqe_temp - fHighPass) * PI / (fHighCut - fHighPass)))
/ ((xishu + 1.0) / (1.0 - xishu) + 1.0);
pSpec[i] *= amp;
}
else
{
pSpec[i] = 0.0;
}
}
if (fLowPass == 0 && fLowPass == fLowCut)
{
pSpec[0] = amp;
}
// 对称
j = L / 2 - 1;
for (i = L / 2; i < L; i++)
{
if (j < 0)
{
break;
}
fqe_temp = i * dalf;
pFreq[i] = fqe_temp;
pSpec[i] = pSpec[j];
j--;
}
return;
}
void DoubleCos_Wavelet(float *pSpec, int nFreqLen, int nK,
int nThickRatio,
int Lr, //中心点
int L, //子波长度
float dalt, //采样间隔 s
float *pTime, float *pValue, float amp)
{
//=========================先在频率域抽稀=============================
int tFreqLen;
float fMaxValue = -1999.0;
tFreqLen = nFreqLen / nThickRatio;
float *pr = new float[tFreqLen]; // 振幅谱
float *pi = new float[tFreqLen]; // 相位谱
float *fr = new float[tFreqLen];
float *fi = new float[tFreqLen];
memset(pr, 0, sizeof(float)*tFreqLen);
int i = 0, j, iCnt = 0;
int Ll, k;
for (i = 0; i < nFreqLen; i += nThickRatio)
{
if (iCnt == tFreqLen)
{
break;
}
pr[iCnt] = pSpec[i];
pi[iCnt] = 0.0;
iCnt++;
}
//float *pr,float *pi,int n,int k,float *fr,float *fi,int l,int il
kkfft(pr, pi, tFreqLen, nK, fr, fi, 1, 0);
Ll = L - Lr - 1;
j = Lr;
for (i = -Lr; i <= 0; i++)
{
k = i + Lr;
pTime[k] = (float)i * dalt * 1000.0;
pValue[k] = amp * fr[j];
if (fMaxValue < pValue[k])
{
fMaxValue = pValue[k];
}
j--;
}
j = L - 1;
for (i = 1; i <= Ll; i++)
{
k = i + Lr;
pTime[k] = (float)i * dalt * 1000.0;
pValue[k] = amp * fr[j];
if (fMaxValue < pValue[k])
{
fMaxValue = pValue[k];
}
j--;
}
if (fMaxValue != 0)
{
for (i = 0; i < L; i++)
{
pValue[i] /= fMaxValue;
}
}
delete[] pr;
delete[] pi;
delete[] fr;
delete[] fi;
return;
}
void Preprocessing(float *Imp_in, unsigned int LenImp, int *start_end) {
for (int i = 0; i < LenImp-1; i++)
{
//出现大于1000的部分记下来
if (Imp_in[i] > Min_Imp) {
start_end[0] = i;
break;
}
}
for (int j = LenImp-1; j >= 0; j--)
{
if (Imp_in[j] > Min_Imp) {
start_end[1] = j;
break;
}
}
}
//计算反射系数
void Compute_Rcs(float *Imp, float* Rcs, int LenImp) {
float **Diff = NULL;
float *LnImp = NULL;
Diff = alloc2float(LenImp - 1, LenImp); //差分矩阵
LnImp = (float*)calloc(LenImp, sizeof(float)); //波阻抗对数值
for (int i = 0; i < LenImp - 1; i++) {
for (int j = 0; j < LenImp; j++) {
if (i == j) { //当矩阵行列下标相等(i=j)时，取值为-0.5
Diff[i][j] = -0.5;
}
else if (i == j - 1) { //当矩阵下标i=j-1时，取值为0.5
Diff[i][j] = 0.5;
}
else { //其他情况，值为0
Diff[i][j] = 0;
}
}
}
for (int i = 0; i < LenImp; i++) {
LnImp[i] = log(Imp[i]); //取对数
}
matrixMultilcation(Diff, LnImp, Rcs, LenImp, LenImp - 1); //矩阵乘法
}
//反射系数褶积子波得到合成记录
void Rcs2Syn(float *Rcs, float *Syn, int LenRcs, int WaveletType, int WaveletLenTemp, int WaveletInc) {
float WaveletAmp = 1.0;
int WaveletLen = WaveletLenTemp /WaveletInc;
float *WaveletValue = NULL;
float *WaveletTime = NULL;
WaveletValue = (float*)calloc(WaveletLen, sizeof(float));
WaveletTime = (float*)calloc(WaveletLen, sizeof(float));
switch (WaveletType) {
case 1:
// type==1 雷克子波
parameterRicker Ricker_para;
Ricker_para.DomainFreq = 30;
Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f,
Ricker_para.DomainFreq, WaveletValue, WaveletTime, WaveletAmp);
break;
case 2:
//type==2 Yu's子波
parameterYus Yus_para;
Yus_para.QFreq = 60.0;
Yus_para.PFreq = 10.0;
Double_Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f, Yus_para.QFreq, Yus_para.PFreq,
WaveletValue, WaveletTime, WaveletAmp);
break;
case 3:
//type==3 双余弦子波
parameterDoubleCos DoubleCos_para;
DoubleCos_para.HighCut = 80;
DoubleCos_para.HighPass = 60;
DoubleCos_para.LowCut = 8;
DoubleCos_para.LowPass = 12;
int nK = getK(WaveletLen);
int SpecLen = (int)pow(2.0, nK - 1);
float FreqInc = 1000.0 / (2.0 * WaveletInc * SpecLen);
int ThickRatio = 8;
int DoubleCosSpecLen = int(SpecLen * 2 * ThickRatio);
float *DoubleCosFreq = NULL;
float *DoubleCosSpec = NULL;
DoubleCosFreq = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCosSpec = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCos_Spec(DoubleCosSpecLen, FreqInc / ThickRatio, DoubleCos_para.LowCut, DoubleCos_para.LowPass,
DoubleCos_para.HighPass, DoubleCos_para.HighCut, DoubleCosFreq, DoubleCosSpec, WaveletAmp);
DoubleCos_Wavelet(DoubleCosSpec, DoubleCosSpecLen, nK, ThickRatio, WaveletLen / 2, WaveletLen, WaveletInc / 1000.0,
WaveletTime, WaveletValue, 1.0);
free(DoubleCosFreq);
free(DoubleCosSpec);
break;
}
float *Signal = NULL;
int lenSignal = WaveletLen + LenRcs - 1;
Signal = (float*)calloc(lenSignal, sizeof(float));
convolve_cwp(LenRcs, 0, Rcs, WaveletLen, 0, WaveletValue, lenSignal, 0, Signal);
memcpy(Syn, Signal + (WaveletLen - 1) / 2, (LenRcs) * sizeof(float));
free(WaveletValue);
free(WaveletTime);
free(Signal);
}
VI main.cpp
#include "ImpToSynProfile.h"
#include"segy.h"
int main(){
const char *filenameInput = "test.segy";
const char *filenameOutput = "Output66.segy";
bhed fileheader;
segy traceheader;
unsigned int nt = 0;
unsigned int sizefileheader = sizeof(fileheader);
unsigned int sizetraceheader = sizeof(traceheader);
unsigned int ncdp = 0;
long long size_file = 0;
long long size_trace = 0;
FILE *fpinput = NULL;
FILE *fpoutput = NULL;
float *dataInput = NULL;
float *dataOutput = NULL;
fpinput = fopen(filenameInput, "rb");
fpoutput = fopen(filenameOutput, "wb");
if (fpinput == NULL) {
printf("can not open %s file\n", filenameInput);
return false;
}
if (fpoutput == NULL) {
printf("can not open %s file\n", filenameOutput);
return false;
}
fread(&fileheader, sizefileheader, 1, fpinput);
nt = fileheader.hns;
_fseeki64(fpinput, 0, SEEK_END);
size_file = _ftelli64(fpinput);
size_trace = nt * sizeof(float) + sizetraceheader;
ncdp = (size_file - (long long)sizefileheader) / size_trace;
_fseeki64(fpinput, sizefileheader, SEEK_SET);
fwrite(&fileheader, sizefileheader, 1, fpoutput);
dataInput = (float*)calloc(nt, sizeof(float));
int WaveletLenTemp = 128;
int WaveletInc = 1;
int WaveletType = 1; //选择子波类型
//遍历每一道
for (int itrace = 0; itrace < ncdp; itrace++)
{
fread(&traceheader, sizetraceheader, 1, fpinput);
fread(dataInput, nt * sizeof(float), 1, fpinput);
//预处理
int *start_end = NULL;
start_end = (int*)calloc(2, sizeof(int));
Preprocessing(dataInput, nt, start_end);
int Start = 0;
int End = 0;
int LenImp = 0;
Start = start_end[0]; //获取有值的部分的起止位置
End = start_end[1];
LenImp = End - Start + 1;
int Boundary = 3; //边界
int LenImp_ex = LenImp + Boundary * 2;
float *Imp = NULL;
float *Rcs = NULL;
float *syn = NULL;
//延展边界2个点
Imp = (float*)calloc(LenImp_ex , sizeof(float));
Rcs = (float*)calloc(LenImp_ex - 1, sizeof(float));
syn = (float*)calloc(LenImp_ex - 1, sizeof(float));
for (int i = Boundary; i < LenImp+Boundary; i++) {
Imp[i] = dataInput[Start-Boundary + i];
}
//上边界
for (int i = 0; i < Boundary; i++) {
Imp[i] = dataInput[Start];
}
//下边界
for (int i = LenImp+Boundary; i < LenImp + Boundary * 2; i++) {
Imp[i] = dataInput[End];
}
//第一步，计算反射系数
Compute_Rcs(Imp,Rcs,LenImp);
//第二步，反射系数褶积子波
Rcs2Syn(Rcs,syn, LenImp_ex - 1, WaveletType, WaveletLenTemp, WaveletInc);
dataOutput = (float*)calloc(nt, sizeof(float));
for (int i = 0; i < LenImp; i++) {
dataOutput[Start + i] = syn[i + Boundary];
// if ((itrace==250)&&(i==250))
// {
// printf("%f", dataOutput[Start+i]);
// }
}
fwrite(&traceheader, sizetraceheader, 1, fpoutput);
fwrite(dataOutput, nt * sizeof(float), 1, fpoutput);
free(start_end);
free(Imp);
free(Rcs);
free(syn);
}
free(dataInput);
free(dataOutput);
fclose(fpinput);
fclose(fpoutput);
return 0;
}
本文共计13469个文字，预计阅读时间需要54分钟。
C++实现地震记录合成+实例：从波阻抗模型中获取与地震对应的反射系数，再将反射系数与子波积分得到合成地震记录。+1+由波阻抗获取反射系数+地震波在介质中传播时，反射系数传递，作用于面。
C/C++实现合成地震记录本实例将从波阻抗模型中获得与之对应的反射系数，再将反射系数与子波褶积得到合成地震记录。
1 由波阻抗获取反射系数地震波在介质中传播时，作用于某个面积上的压力与单位时间内垂直通过此面积的质点流量(即面积乘质点振动速度）之比，具有阻力的含义，称为波阻抗，其数值等于介质密度\(\rho\)与波速\(v\)的乘积，即\(Z=\rho v\)。波阻抗差异形成地震反射，用反射系数表征上下界面差异，公式为
\[{r} = \frac{{{Z_2} - {Z_1}}}{{{Z_2} + {Z_1}}} = \frac{{{\rho _2}{v_2} - {\rho _1}{v_1}}}{{{\rho _2}{v_2} + {\rho _1}{v_1}}}\\Z_2=Z_1\frac{{1 + r}}{{1 - r}}\\ Z_{2}=Z_{1} \exp \left[\ln \left(\frac{1+r_{}}{1-r_{}}\right)\right] \]在实际情况下，反射系数\(r\)一般远远小于 1，即 \(r \ll 1\) 。则有如下关系：
\[\ln \left(\frac{1+r_{}}{1-r_{}}\right) \approx 2 r_{} \]上式可化为：
\[Z_{2} \approx Z_{1} e^{2 r} \\ r_1=\frac{1}{2}(\ln Z_2-\ln Z_1) \]于是，通过对上下两层波阻抗的对数值作差，就能由此得到对应的反射系数。
将上式用矩阵的形式表示为：
\[{\left[ {\begin{array}{*{20}{c}} { - \frac{1}{2}}&{\frac{1}{2}}&{}&{}&{}&{}\\ {}&{ - \frac{1}{2}}&{\frac{1}{2}}&{}&{}&{}\\ {}&{}& \ddots & \ddots &{}&{}\\ {}&{}&{}&{ - \frac{1}{2}}&{\frac{1}{2}}&{}\\ {}&{}&{}&{}&{ - \frac{1}{2}}&{\frac{1}{2}} \end{array}} \right]_{(n - 1) \times n}}\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\ln Z_1}}\\ {{\ln Z_2}} \end{array}}\\ \vdots \\ {{\ln Z_{n-1}}} \end{array}}\\ {{\ln Z_{n}}} \end{array}} \right]_{n\times 1} = \left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{r_1}}\\ {{r_1}} \end{array}}\\ \vdots \\ {{r_{n-2}}} \end{array}}\\ {{r_{n-1}}} \end{array}} \right]_{(n-1)\times 1} \]容易看出整个实现过程，先是对波阻抗序列求对数，之后左乘差分矩阵，便能够得到反射系数矩阵（长度=\(n-1\)）
代码如下：
/*
输入：
波阻抗 float* Imp
波阻抗序列长度 int nt
输出：
反射系数 float* Rcs
*/
void Compute_Rcs(float *Imp, float* Rcs, int LenImp) {
float **Diff = NULL;
float *LnImp = NULL;
Diff = alloc2float(LenImp - 1, LenImp); //差分矩阵
LnImp = (float*)calloc(LenImp, sizeof(float)); //波阻抗对数值
for (int i = 0; i < LenImp - 1; i++) {
for (int j = 0; j < LenImp; j++) {
if (i == j) { //当矩阵行列下标相等(i=j)时，取值为-0.5
Diff[i][j] = -0.5;
}
else if (i == j - 1) { //当矩阵下标i=j-1时，取值为0.5
Diff[i][j] = 0.5;
}
else { //其他情况，值为0
Diff[i][j] = 0;
}
}
}
for (int i = 0; i < LenImp; i++) {
LnImp[i] = log(Imp[i]); //取对数
}
matrixMultilcation(Diff, LnImp, Rcs, LenImp, LenImp - 1); //矩阵乘法
}
2 褶积模型
地震道可以用褶积模型描述：
\[S(t) =w(t) *r(t) \]其中\(S(t)\)为地震道，\(w(t)\)为地震子波，\(r(t)\)为反射系数。震源在地表激发的时候震源在地表激发的时刻，地震子波被认为是一个尖脉冲信号。如果这个脉冲在传播过程中能量和波形均保持不变，那么地震记录就是反射系数剖面。而实际上子波在传播过程中，由于受球面扩散和吸收衰减的影响，并非严格的尖脉冲信号。
在本次实例中，使用的子波有以下三种：
/*
函数功能：
-由反射系数褶积地震子波，得到合成地震记录
输入：
-反射系数序列 float *Rcs
-反射系数长度 int LenRcs
-子波类型 int WaveletType
-子波长度 int WaveletLenTemp
-子波采样间隔 int WaveletInc
输出：
-合成地震记录 float *Syn
*/
void Rcs2Syn(float *Rcs, float *Syn, int LenRcs, int WaveletType, int WaveletLenTemp, int WaveletInc) {
float WaveletAmp = 1.0;
int WaveletLen = WaveletLenTemp /WaveletInc;
float *WaveletValue = NULL;
float *WaveletTime = NULL;
WaveletValue = (float*)calloc(WaveletLen, sizeof(float));
WaveletTime = (float*)calloc(WaveletLen, sizeof(float));
switch (WaveletType) {
case 1:
// type==1 雷克子波
parameterRicker Ricker_para;
Ricker_para.DomainFreq = 30;
Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f,
Ricker_para.DomainFreq, WaveletValue, WaveletTime, WaveletAmp);
break;
case 2:
//type==2 Yu's子波
parameterYus Yus_para;
Yus_para.QFreq = 60.0;
Yus_para.PFreq = 10.0;
Double_Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f, Yus_para.QFreq, Yus_para.PFreq,
WaveletValue, WaveletTime, WaveletAmp);
break;
case 3:
//type==3 双余弦子波
parameterDoubleCos DoubleCos_para;
DoubleCos_para.HighCut = 80;
DoubleCos_para.HighPass = 60;
DoubleCos_para.LowCut = 8;
DoubleCos_para.LowPass = 12;
int nK = getK(WaveletLen);
int SpecLen = (int)pow(2.0, nK - 1);
float FreqInc = 1000.0 / (2.0 * WaveletInc * SpecLen);
int ThickRatio = 8;
int DoubleCosSpecLen = int(SpecLen * 2 * ThickRatio);
float *DoubleCosFreq = NULL;
float *DoubleCosSpec = NULL;
DoubleCosFreq = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCosSpec = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCos_Spec(DoubleCosSpecLen, FreqInc / ThickRatio, DoubleCos_para.LowCut, DoubleCos_para.LowPass,
DoubleCos_para.HighPass, DoubleCos_para.HighCut, DoubleCosFreq, DoubleCosSpec, WaveletAmp);
DoubleCos_Wavelet(DoubleCosSpec, DoubleCosSpecLen, nK, ThickRatio, WaveletLen / 2, WaveletLen, WaveletInc / 1000.0,
WaveletTime, WaveletValue, 1.0);
free(DoubleCosFreq);
free(DoubleCosSpec);
break;
}
float *Signal = NULL;
int lenSignal = WaveletLen + LenRcs - 1;
Signal = (float*)calloc(lenSignal, sizeof(float));
convolve_cwp(LenRcs, 0, Rcs, WaveletLen, 0, WaveletValue, lenSignal, 0, Signal); //卷积
memcpy(Syn, Signal + (WaveletLen - 1) / 2, (LenRcs) * sizeof(float)); //赋值
free(WaveletValue);
free(WaveletTime);
free(Signal);
}
3 主函数
输入的波阻抗数据将被转换为合成记录：
int main(){
const char *filenameInput = "Imp.segy"; //输入segy文件
const char *filenameOutput = "Syn.segy"; //输出segy文件
bhed fileheader;
segy traceheader;
unsigned int nt = 0;
unsigned int sizefileheader = sizeof(fileheader);
unsigned int sizetraceheader = sizeof(traceheader);
unsigned int ncdp = 0;
long long size_file = 0;
long long size_trace = 0;
FILE *fpinput = NULL;
FILE *fpoutput = NULL;
float *dataInput = NULL;
float *dataOutput = NULL;
fpinput = fopen(filenameInput, "rb");
fpoutput = fopen(filenameOutput, "wb");
if (fpinput == NULL) {
printf("can not open %s file\n", filenameInput);
return false;
}
if (fpoutput == NULL) {
printf("can not open %s file\n", filenameOutput);
return false;
}
fread(&fileheader, sizefileheader, 1, fpinput);
nt = fileheader.hns;
_fseeki64(fpinput, 0, SEEK_END);
size_file = _ftelli64(fpinput);
size_trace = nt * sizeof(float) + sizetraceheader;
ncdp = (size_file - (long long)sizefileheader) / size_trace;
_fseeki64(fpinput, sizefileheader, SEEK_SET);
fwrite(&fileheader, sizefileheader, 1, fpoutput);
dataInput = (float*)calloc(nt, sizeof(float));
/****************以上为segy文件的读入*****************/
int WaveletLenTemp = 128; //子波长度
int WaveletInc = 1; //子波间隔
int WaveletType = 1; //选择子波类型，此处为Ricker
//遍历每一道
for (int itrace = 0; itrace < ncdp; itrace++)
{
fread(&traceheader, sizetraceheader, 1, fpinput);
fread(dataInput, nt * sizeof(float), 1, fpinput);
//预处理
int *start_end = NULL;
start_end = (int*)calloc(2, sizeof(int));
Preprocessing(dataInput, nt, start_end);//Preprocessing的作用是屏蔽0值和明显的异常值，读取有值的部分
int Start = 0;
int End = 0;
int LenImp = 0;
Start = start_end[0]; //获取有值部分的开始和结尾
End = start_end[1];
LenImp = End - Start + 1;
int Boundary = 3; //为了防止边界出现假轴，将波阻抗稍微延长出一部分
int LenImp_ex = LenImp + Boundary * 2;
float *Imp = NULL;
float *Rcs = NULL;
float *syn = NULL;
Imp = (float*)calloc(LenImp_ex , sizeof(float));
Rcs = (float*)calloc(LenImp_ex - 1, sizeof(float));
syn = (float*)calloc(LenImp_ex - 1, sizeof(float));
for (int i = Boundary; i < LenImp+Boundary; i++) {
Imp[i] = dataInput[Start-Boundary + i];
}
//上边界
for (int i = 0; i < Boundary; i++) {
Imp[i] = dataInput[Start];
}
//下边界
for (int i = LenImp+Boundary; i < LenImp + Boundary * 2; i++) {
Imp[i] = dataInput[End];
}
//第一步，计算反射系数
Compute_Rcs(Imp,Rcs,LenImp);
//第二步，反射系数褶积子波，得到合成记录
Rcs2Syn(Rcs,syn, LenImp_ex - 1, WaveletType, WaveletLenTemp, WaveletInc);
dataOutput = (float*)calloc(nt, sizeof(float));
for (int i = 0; i < LenImp; i++) {
dataOutput[Start + i] = syn[i + Boundary]; //传入合成记录
}
//写入segy
fwrite(&traceheader, sizetraceheader, 1, fpoutput);
fwrite(dataOutput, nt * sizeof(float), 1, fpoutput);
free(start_end);
free(Imp);
free(Rcs);
free(syn);
}
free(dataInput);
free(dataOutput);
fclose(fpinput);
fclose(fpoutput);
return 0;
}
完整代码
I fundamental.h
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include<iostream>
#include<vector>
#include<algorithm>
#include"complex.h"
#include"alloc.h"
void matrixMultilcation(float **data1, float *data2, float *dataOut, int c, int m); //矩阵乘法
void convolve_cwp(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z); //褶积
int getK(int N); //求2的K次方
void kkfft(float *pr, float *pi, int n, int k, float *fr, float *fi, int l, int il);//傅氏变换
II fundamental.cpp
#include "fundamental.h"
// 矩阵乘法
// data1 m行c列
//data2 c行1列
void matrixMultilcation(float **data1, float *data2, float *dataOut, int c, int m) {
float sumTemp = 0;
for (int i = 0; i < m; i++) {
sumTemp = 0;
for (int k = 0; k < c; k++) {
sumTemp += data1[i][k] * data2[k];
}
dataOut[i] = sumTemp;
}
}
// 满足2^k>=N的最小k值,返回k值
int getK(int N)
{
int k, temp = 1;
for (k = 1; k < N; k++)
{
temp *= 2;
if (temp >= N)
{
break;
}
}
return k;
}
void kkfft(float *pr, float *pi, int n, int k, float *fr, float *fi, int l, int il)
{
int it, m, is, i, j, nv, l0;
float p, q, s, vr, vi, poddr, poddi;
for (it = 0; it <= n - 1; it++)
{
m = it;
is = 0;
for (i = 0; i <= k - 1; i++)
{
j = m / 2;
is = 2 * is + (m - 2 * j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306 / (1.0 * n);
pr[1] = cos(p);
pi[1] = -sin(p);
if (l != 0)
{
pi[1] = -pi[1];
}
for (i = 2; i <= n - 1; i++)
{
p = pr[i - 1] * pr[1];
q = pi[i - 1] * pi[1];
s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
pr[i] = p - q;
pi[i] = s - p - q;
}
for (it = 0; it <= n - 2; it = it + 2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr + fr[it + 1];
fi[it] = vi + fi[it + 1];
fr[it + 1] = vr - fr[it + 1];
fi[it + 1] = vi - fi[it + 1];
}
m = n / 2;
nv = 2;
for (l0 = k - 2; l0 >= 0; l0--)
{
m = m / 2;
nv = 2 * nv;
for (it = 0; it <= (m - 1)*nv; it = it + nv)
{
for (j = 0; j <= (nv / 2) - 1; j++)
{
p = pr[m * j] * fr[it + j + nv / 2];
q = pi[m * j] * fi[it + j + nv / 2];
s = pr[m * j] + pi[m * j];
s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
poddr = p - q;
poddi = s - p - q;
fr[it + j + nv / 2] = fr[it + j] - poddr;
fi[it + j + nv / 2] = fi[it + j] - poddi;
fr[it + j] = fr[it + j] + poddr;
fi[it + j] = fi[it + j] + poddi;
}
}
}
if (l != 0)
{
for (i = 0; i <= n - 1; i++)
{
fr[i] = fr[i] / (1.0 * n);
fi[i] = fi[i] / (1.0 * n);
}
if (il != 0)
{
for (i = 0; i <= n - 1; i++)
{
pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]);
if (fabs(fr[i]) < 0.000001 * fabs(fi[i]))
{
if ((fi[i] * fr[i]) > 0)
{
pi[i] = 90.0;
}
else
{
pi[i] = -90.0;
}
}
else
{
pi[i] = atan(fi[i] / fr[i]) * 360.0 / 6.283185306;
}
}
}
}
return;
}
III convolution.cpp
/******************************************************************************
Input:
lx length of x array
ifx sample index of first x
x array[lx] to be convolved with y
ly length of y array
ify sample index of first y
y array[ly] with which x is to be convolved
lz length of z array
ifz sample index of first z
Output:
z array[lz] containing x convolved with y
******************************************************************************/
#include"fundamental.h"
/* prototype of function used internally */
static void convolve_cwp_s(int, int, float*, int, int, float*, int, int, float*);
void convolve_cwp(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z)
#ifdef SIMPLE_CONV
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1, i, j, jlow, jhigh;
float sum;
x -= ifx; y -= ify; z -= ifz;
for (i = ifz; i <= ilz; ++i) {
jlow = i - ily; if (jlow < ifx) jlow = ifx;
jhigh = i - ify; if (jhigh > ilx) jhigh = ilx;
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
}
#else
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1,
i, j, ilow, ihigh, jlow, jhigh;
float sa, sb, xa, xb, ya, yb, *t;
/* if x is longer than y, swap x and y */
if (lx > ly) {
i = ifx; ifx = ify; ify = i;
i = ilx; ilx = ily; ily = i;
i = lx; lx = ly; ly = i;
t = x; x = y; y = t;
}
/* handle short x with special code */
if (lx >= 1 && lx <= 30) {
convolve_cwp_s(lx, ifx, x, ly, ify, y, lz, ifz, z);
return;
}
/* adjust pointers for indices of first samples */
x -= ifx;
y -= ify;
z -= ifz;
/* OFF LEFT: i < ify+ifx */
/* zero output for all i */
ilow = ifz;
ihigh = ify + ifx - 1; if (ihigh > ilz) ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
/* ROLLING ON: ify+ifx <= i < ify+ilx */
/* if necessary, do one i so that number of j in overlap is odd */
if (i < ify + ilx && i <= ilz) {
jlow = ifx;
jhigh = i - ify;
if ((jhigh - jlow) % 2) {
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
}
/* loop over pairs of i and j */
ilow = i;
ihigh = ilx + ify - 1; if (ihigh > ilz) ihigh = ilz;
jlow = ifx;
jhigh = ilow - ify;
for (i = ilow; i < ihigh; i += 2, jhigh += 2) {
sa = sb = 0.0;
xb = x[jhigh + 1];
yb = 0.0;
for (j = jhigh; j >= jlow; j -= 2) {
sa += xb * yb;
ya = y[i - j];
sb += xb * ya;
xa = x[j];
sa += xa * ya;
yb = y[i + 1 - j];
sb += xa * yb;
xb = x[j - 1];
}
z[i] = sa;
z[i + 1] = sb;
}
/* if number of i is odd */
if (i == ihigh) {
jlow = ifx;
jhigh = i - ify;
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
/* MIDDLE: ify+ilx <= i <= ily+ifx */
/* determine limits for i and j */
ilow = i;
ihigh = ily + ifx; if (ihigh > ilz) ihigh = ilz;
jlow = ifx;
jhigh = ilx;
/* if number of j is even, do j in pairs with no leftover */
if ((jhigh - jlow) % 2) {
for (i = ilow; i < ihigh; i += 2) {
sa = sb = 0.0;
yb = y[i + 1 - jlow];
xa = x[jlow];
for (j = jlow; j < jhigh; j += 2) {
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
xa = x[j + 2];
}
z[i] = sa;
z[i + 1] = sb;
}
/* else, number of j is odd, so do j in pairs with leftover */
}
else {
for (i = ilow; i < ihigh; i += 2) {
sa = sb = 0.0;
yb = y[i + 1 - jlow];
xa = x[jlow];
for (j = jlow; j < jhigh; j += 2) {
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
xa = x[j + 2];
}
z[i] = sa + x[jhigh] * y[i - jhigh];
z[i + 1] = sb + x[jhigh] * y[i + 1 - jhigh];
}
}
/* if number of i is odd */
if (i == ihigh) {
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
/* ROLLING OFF: ily+ifx < i <= ily+ilx */
/* if necessary, do one i so that number of j in overlap is even */
if (i <= ily + ilx && i <= ilz) {
jlow = i - ily;
jhigh = ilx;
if (!((jhigh - jlow) % 2)) {
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
}
/* number of j is now even, so loop over both i and j in pairs */
ilow = i;
ihigh = ily + ilx; if (ihigh > ilz) ihigh = ilz;
jlow = ilow - ily;
jhigh = ilx - 2; /* Dave's new patch */
for (i = ilow; i < ihigh; i += 2, jlow += 2) {
sa = sb = 0.0;
xa = x[jlow];
yb = 0.0;
for (j = jlow; j < jhigh; j += 2) {
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
xa = x[j + 2];
}
sb += xa * yb;
ya = y[i - j];
sa += xa * ya;
xb = x[j + 1];
sb += xb * ya;
yb = y[i - 1 - j];
sa += xb * yb;
z[i] = sa;
z[i + 1] = sb;
}
/* if number of i is odd */
if (i == ihigh) {
jlow = i - ily;
jhigh = ilx;
sa = 0.0;
for (j = jlow; j <= jhigh; ++j)
sa += x[j] * y[i - j];
z[i++] = sa;
}
/* OFF RIGHT: ily+ilx < i */
/* zero output for all i */
ilow = i;
ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
}
/* internal function optimized for short x */
static void convolve_cwp_s(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z)
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1,
i, j, ilow, ihigh, jlow, jhigh;
float x0, x1, x2, x3, x4, x5, x6, x7, x8, x9,
x10, x11, x12, x13, x14, x15, x16, x17, x18, x19,
x20, x21, x22, x23, x24, x25, x26, x27, x28, x29,
ya, yb, z0, z1, sum;
x -= ifx;
y -= ify;
z -= ifz;
/* OFF LEFT: i < ifx+ify */
ilow = ifz;
ihigh = ify + ifx - 1; if (ihigh > ilz) ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
/* ROLLING ON: ify+ifx <= i < ify+ilx */
ilow = ify + ifx; if (ilow < ifz) ilow = ifz;
ihigh = ify + ilx - 1; if (ihigh > ilz) ihigh = ilz;
jlow = ifx;
jhigh = ilow - ify;
for (i = ilow; i <= ihigh; ++i, ++jhigh) {
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
/* MIDDLE: ify+ilx <= i <= ily+ifx */
ilow = ify + ilx; if (ilow < ifz) ilow = ifz;
ihigh = ily + ifx; if (ihigh > ilz) ihigh = ilz;
if (lx == 1) {
x0 = x[ifx];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 2) {
x0 = x[ifx];
x1 = x[ifx + 1];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 3) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 4) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 5) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 6) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 7) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 8) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 9) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 10) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 11) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 12) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 13) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 14) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 15) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 16) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 17) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 18) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 19) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 20) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 21) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 22) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 23) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 24) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 25) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 26) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 27) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 28) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
x27 = x[ifx + 27];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb; z1 += x27 * yb;
ya = y[i - ifx - 27]; z0 += x27 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 29) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
x27 = x[ifx + 27];
x28 = x[ifx + 28];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb; z1 += x27 * yb;
ya = y[i - ifx - 27]; z0 += x27 * ya; z1 += x28 * ya;
yb = y[i - ifx - 28]; z0 += x28 * yb;
z[i + 1] = z1;
z[i] = z0;
}
}
else if (lx == 30) {
x0 = x[ifx];
x1 = x[ifx + 1];
x2 = x[ifx + 2];
x3 = x[ifx + 3];
x4 = x[ifx + 4];
x5 = x[ifx + 5];
x6 = x[ifx + 6];
x7 = x[ifx + 7];
x8 = x[ifx + 8];
x9 = x[ifx + 9];
x10 = x[ifx + 10];
x11 = x[ifx + 11];
x12 = x[ifx + 12];
x13 = x[ifx + 13];
x14 = x[ifx + 14];
x15 = x[ifx + 15];
x16 = x[ifx + 16];
x17 = x[ifx + 17];
x18 = x[ifx + 18];
x19 = x[ifx + 19];
x20 = x[ifx + 20];
x21 = x[ifx + 21];
x22 = x[ifx + 22];
x23 = x[ifx + 23];
x24 = x[ifx + 24];
x25 = x[ifx + 25];
x26 = x[ifx + 26];
x27 = x[ifx + 27];
x28 = x[ifx + 28];
x29 = x[ifx + 29];
for (i = ilow; i <= ihigh - 1; i += 2) {
ya = y[i + 1 - ifx]; z1 = x0 * ya;
yb = y[i - ifx]; z0 = x0 * yb; z1 += x1 * yb;
ya = y[i - ifx - 1]; z0 += x1 * ya; z1 += x2 * ya;
yb = y[i - ifx - 2]; z0 += x2 * yb; z1 += x3 * yb;
ya = y[i - ifx - 3]; z0 += x3 * ya; z1 += x4 * ya;
yb = y[i - ifx - 4]; z0 += x4 * yb; z1 += x5 * yb;
ya = y[i - ifx - 5]; z0 += x5 * ya; z1 += x6 * ya;
yb = y[i - ifx - 6]; z0 += x6 * yb; z1 += x7 * yb;
ya = y[i - ifx - 7]; z0 += x7 * ya; z1 += x8 * ya;
yb = y[i - ifx - 8]; z0 += x8 * yb; z1 += x9 * yb;
ya = y[i - ifx - 9]; z0 += x9 * ya; z1 += x10 * ya;
yb = y[i - ifx - 10]; z0 += x10 * yb; z1 += x11 * yb;
ya = y[i - ifx - 11]; z0 += x11 * ya; z1 += x12 * ya;
yb = y[i - ifx - 12]; z0 += x12 * yb; z1 += x13 * yb;
ya = y[i - ifx - 13]; z0 += x13 * ya; z1 += x14 * ya;
yb = y[i - ifx - 14]; z0 += x14 * yb; z1 += x15 * yb;
ya = y[i - ifx - 15]; z0 += x15 * ya; z1 += x16 * ya;
yb = y[i - ifx - 16]; z0 += x16 * yb; z1 += x17 * yb;
ya = y[i - ifx - 17]; z0 += x17 * ya; z1 += x18 * ya;
yb = y[i - ifx - 18]; z0 += x18 * yb; z1 += x19 * yb;
ya = y[i - ifx - 19]; z0 += x19 * ya; z1 += x20 * ya;
yb = y[i - ifx - 20]; z0 += x20 * yb; z1 += x21 * yb;
ya = y[i - ifx - 21]; z0 += x21 * ya; z1 += x22 * ya;
yb = y[i - ifx - 22]; z0 += x22 * yb; z1 += x23 * yb;
ya = y[i - ifx - 23]; z0 += x23 * ya; z1 += x24 * ya;
yb = y[i - ifx - 24]; z0 += x24 * yb; z1 += x25 * yb;
ya = y[i - ifx - 25]; z0 += x25 * ya; z1 += x26 * ya;
yb = y[i - ifx - 26]; z0 += x26 * yb; z1 += x27 * yb;
ya = y[i - ifx - 27]; z0 += x27 * ya; z1 += x28 * ya;
yb = y[i - ifx - 28]; z0 += x28 * yb; z1 += x29 * yb;
ya = y[i - ifx - 29]; z0 += x29 * ya;
z[i + 1] = z1;
z[i] = z0;
}
}
if (ihigh >= ilow && (ihigh - ilow) % 2 == 0) {
ilow = ihigh;
jlow = ifx;
jhigh = ilx;
for (i = ilow; i <= ihigh; ++i) {
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
}
/* ROLLING OFF: ily+ifx < i <= ily+ilx */
ilow = ily + ifx + 1; if (ilow < ifz) ilow = ifz;
ihigh = ily + ilx; if (ihigh > ilz) ihigh = ilz;
jlow = ilow - ily;
jhigh = ilx;
for (i = ilow; i <= ihigh; ++i, ++jlow) {
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
/* OFF RIGHT: ily+ilx < i */
ilow = ily + ilx + 1; if (ilow < ifz) ilow = ifz;
ihigh = ilz;
for (i = ilow; i <= ihigh; ++i)
z[i] = 0.0;
}
#endif
#ifdef TEST
static void convolve_cwp1(int lx, int ifx, float *x,
int ly, int ify, float *y,
int lz, int ifz, float *z)
{
int ilx = ifx + lx - 1, ily = ify + ly - 1, ilz = ifz + lz - 1, i, j, jlow, jhigh;
float sum;
x -= ifx; y -= ify; z -= ifz;
for (i = ifz; i <= ilz; ++i) {
jlow = i - ily; if (jlow < ifx) jlow = ifx;
jhigh = i - ify; if (jhigh > ilx) jhigh = ilx;
for (j = jlow, sum = 0.0; j <= jhigh; ++j)
sum += x[j] * y[i - j];
z[i] = sum;
}
}
#include "cwp.h"
main()
{
int lx, ly, lz, ifx, ify, ifz, i;
float *x, *y, *z, *z1;
/* loop over tests */
while (1) {
/* make parameters */
lx = 1 + 100 * franuni();
ly = 1 + 100 * franuni();
lz = 1 + 100 * franuni();
ifx = -100 + 200 * franuni();
ify = -100 + 200 * franuni();
ifz = -100 + 200 * franuni();
/* allocate space */
x = alloc1float(lx);
y = alloc1float(ly);
z = alloc1float(lz);
z1 = alloc1float(lz);
/* fill x and y arrays */
for (i = 0; i < lx; i++)
x[i] = i + 1;
for (i = 0; i < ly; i++)
y[i] = i + 1;
/* convolve and check */
convolve_cwp1(lx, ifx, x, ly, ify, y, lz, ifz, z1);
convolve_cwp(lx, ifx, x, ly, ify, y, lz, ifz, z);
for (i = 0; i < lz; ++i)
if (z[i] != z1[i]) break;
if (i < lz) {
printf("%10d %10d %10d %10d %10d %10d\n",
lx, ifx, ly, ify, lz, ifz);
printf("z1[%d]=%g != z[%d]=%g\n",
i + ifz, z1[i], i + ifz, z[i]);
pp1d(stdout, "simple", lz, ifz, z1);
pp1d(stdout, "optimized", lz, ifz, z);
exit(-1);
}
/* free space */
free1float(x);
free1float(y);
free1float(z);
free1float(z1);
}
}
#endif
IV ImpToSyn.h
#include"fundamental.h"
using namespace std;
//雷克子波
void Ricker(int Lr, //中心点
int L, //子波长度
float dalt, //采样间隔，单位：s
float fp, //主频
float *w, //子波数组
float *time,//横轴数组
float amp); //最大振幅
void Double_Ricker(int Lr, int L, float dalt, float fq, float fp, float *w, float *time, float amp);//Yu'S
void DoubleCos_Spec(int L, float dalf, float fLowCut, float fLowPass, float fHighPass, float fHighCut, float *pFreq, float *pSpec, float amp);//双余弦
void DoubleCos_Wavelet(float *pSpec, int nFreqLen, int nK,int nThickRatio,int Lr,int L,float dalt,float *pTime, float *pValue, float amp);
void Compute_Rcs(float *Imp, float* Rcs,int LenImp); //计算反射系数
void Rcs2Syn(float *Rcs, float *Syn, int LenRcs, int WaveletType, int WaveletLenTemp,int WaveletInc); //反射系数褶积子波得到合成记录
void Preprocessing(float *Imp_in, unsigned int LenImp, int *start_end); //对输入地震道进行预处理
typedef struct parameterwavelet {
int type; // 子波类型
int WaveletLen;
int WaveletInc;
}parameterwavelet;
typedef struct parameterRicker { // 雷克子波参数
int DomainFreq;
}parameterRicker;
typedef struct parameterYus { // Yu's子波参数
float QFreq;
float PFreq;
}parameterYus;
typedef struct parameterDoubleCos { // DoubleCos子波参数
float LowCut;
float LowPass;
float HighPass;
float HighCut;
}parameterDoubleCos;
V ImpToSyn.cpp
#include"ImpToSynProfile.h"
#define Min_Imp 1000.0
/*雷克子波*/
void Ricker(int Lr, //中心点
int L, //子波长度
float dalt, //采样间隔，单位：s
float fp, //主频
float *w, //子波数组
float *time,//横轴数组
float amp) //最大振幅
{
double temp;
double pi;
int Ll;
int k, i;
Ll = L - Lr - 1;
pi = 4.0 * atan(1.0);
for (i = -Lr; i <= Ll; i++)
{
k = i + Lr;
temp = pi * fp * (dalt * i);
temp = temp * temp;
time[k] = (float)i * dalt * 1000.0;
w[k] = amp * (1.0 - 2.0 * temp) * exp(-temp);
}
return;
}
/*Yu's子波*/
void Double_Ricker(int Lr, int L, float dalt, float fq, float fp, float *w, float *time, float amp) //与Ricker相同，有Q,P两个频率
{
double temp1, temp2;
double pi;
int Ll;
int k, i;
Ll = L - Lr - 1;
pi = 4.0 * atan(1.0);
for (i = -Lr; i <= Ll; i++)
{
k = i + Lr;
temp1 = pi * fq * (dalt * i);
temp2 = pi * fp * (dalt * i);
temp1 = temp1 * temp1;
temp2 = temp2 * temp2;
time[k] = (float)i * dalt * 1000.0;
w[k] = amp * (1.0 / (fq - fp)) * (fq * exp(-temp1) - fp * exp(-temp2));
}
return;
}
//计算双余弦子波频谱
void DoubleCos_Spec(int L, float dalf, float fLowCut, float fLowPass, float fHighPass, float fHighCut, float *pFreq, float *pSpec, float amp)
{
//------ bak ------
int i, j;
float fqe_temp;
float PI;
float xishu = 0;
if (L % 2 > 0)
{
L--;
}
PI = 4.0 * atan(1.0);
for (i = 0; i < L / 2; i++)
{
fqe_temp = i * dalf;
pFreq[i] = fqe_temp;
if (fqe_temp <= fLowCut)
{
pSpec[i] = 0.0;
}
else if (fLowCut < fLowPass && fqe_temp > fLowCut && fqe_temp < fLowPass)
{
pSpec[i] = 1.0 * ((xishu + 1.0) / (1.0 - xishu)
+ cos((fqe_temp - fLowCut) * PI / (fLowPass - fLowCut) + PI))
/ ((xishu + 1.0) / (1.0 - xishu) + 1.0);
pSpec[i] *= amp;
}
else if (fqe_temp >= fLowPass && fqe_temp <= fHighPass)
{
pSpec[i] = 1.0;
pSpec[i] *= amp;
}
else if (fHighPass < fHighCut && fqe_temp > fHighPass && fqe_temp < fHighCut)
{
pSpec[i] = 1.0 * ((xishu + 1.0) / (1.0 - xishu)
+ cos((fqe_temp - fHighPass) * PI / (fHighCut - fHighPass)))
/ ((xishu + 1.0) / (1.0 - xishu) + 1.0);
pSpec[i] *= amp;
}
else
{
pSpec[i] = 0.0;
}
}
if (fLowPass == 0 && fLowPass == fLowCut)
{
pSpec[0] = amp;
}
// 对称
j = L / 2 - 1;
for (i = L / 2; i < L; i++)
{
if (j < 0)
{
break;
}
fqe_temp = i * dalf;
pFreq[i] = fqe_temp;
pSpec[i] = pSpec[j];
j--;
}
return;
}
void DoubleCos_Wavelet(float *pSpec, int nFreqLen, int nK,
int nThickRatio,
int Lr, //中心点
int L, //子波长度
float dalt, //采样间隔 s
float *pTime, float *pValue, float amp)
{
//=========================先在频率域抽稀=============================
int tFreqLen;
float fMaxValue = -1999.0;
tFreqLen = nFreqLen / nThickRatio;
float *pr = new float[tFreqLen]; // 振幅谱
float *pi = new float[tFreqLen]; // 相位谱
float *fr = new float[tFreqLen];
float *fi = new float[tFreqLen];
memset(pr, 0, sizeof(float)*tFreqLen);
int i = 0, j, iCnt = 0;
int Ll, k;
for (i = 0; i < nFreqLen; i += nThickRatio)
{
if (iCnt == tFreqLen)
{
break;
}
pr[iCnt] = pSpec[i];
pi[iCnt] = 0.0;
iCnt++;
}
//float *pr,float *pi,int n,int k,float *fr,float *fi,int l,int il
kkfft(pr, pi, tFreqLen, nK, fr, fi, 1, 0);
Ll = L - Lr - 1;
j = Lr;
for (i = -Lr; i <= 0; i++)
{
k = i + Lr;
pTime[k] = (float)i * dalt * 1000.0;
pValue[k] = amp * fr[j];
if (fMaxValue < pValue[k])
{
fMaxValue = pValue[k];
}
j--;
}
j = L - 1;
for (i = 1; i <= Ll; i++)
{
k = i + Lr;
pTime[k] = (float)i * dalt * 1000.0;
pValue[k] = amp * fr[j];
if (fMaxValue < pValue[k])
{
fMaxValue = pValue[k];
}
j--;
}
if (fMaxValue != 0)
{
for (i = 0; i < L; i++)
{
pValue[i] /= fMaxValue;
}
}
delete[] pr;
delete[] pi;
delete[] fr;
delete[] fi;
return;
}
void Preprocessing(float *Imp_in, unsigned int LenImp, int *start_end) {
for (int i = 0; i < LenImp-1; i++)
{
//出现大于1000的部分记下来
if (Imp_in[i] > Min_Imp) {
start_end[0] = i;
break;
}
}
for (int j = LenImp-1; j >= 0; j--)
{
if (Imp_in[j] > Min_Imp) {
start_end[1] = j;
break;
}
}
}
//计算反射系数
void Compute_Rcs(float *Imp, float* Rcs, int LenImp) {
float **Diff = NULL;
float *LnImp = NULL;
Diff = alloc2float(LenImp - 1, LenImp); //差分矩阵
LnImp = (float*)calloc(LenImp, sizeof(float)); //波阻抗对数值
for (int i = 0; i < LenImp - 1; i++) {
for (int j = 0; j < LenImp; j++) {
if (i == j) { //当矩阵行列下标相等(i=j)时，取值为-0.5
Diff[i][j] = -0.5;
}
else if (i == j - 1) { //当矩阵下标i=j-1时，取值为0.5
Diff[i][j] = 0.5;
}
else { //其他情况，值为0
Diff[i][j] = 0;
}
}
}
for (int i = 0; i < LenImp; i++) {
LnImp[i] = log(Imp[i]); //取对数
}
matrixMultilcation(Diff, LnImp, Rcs, LenImp, LenImp - 1); //矩阵乘法
}
//反射系数褶积子波得到合成记录
void Rcs2Syn(float *Rcs, float *Syn, int LenRcs, int WaveletType, int WaveletLenTemp, int WaveletInc) {
float WaveletAmp = 1.0;
int WaveletLen = WaveletLenTemp /WaveletInc;
float *WaveletValue = NULL;
float *WaveletTime = NULL;
WaveletValue = (float*)calloc(WaveletLen, sizeof(float));
WaveletTime = (float*)calloc(WaveletLen, sizeof(float));
switch (WaveletType) {
case 1:
// type==1 雷克子波
parameterRicker Ricker_para;
Ricker_para.DomainFreq = 30;
Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f,
Ricker_para.DomainFreq, WaveletValue, WaveletTime, WaveletAmp);
break;
case 2:
//type==2 Yu's子波
parameterYus Yus_para;
Yus_para.QFreq = 60.0;
Yus_para.PFreq = 10.0;
Double_Ricker(WaveletLen / 2, WaveletLen, WaveletInc / 1000.f, Yus_para.QFreq, Yus_para.PFreq,
WaveletValue, WaveletTime, WaveletAmp);
break;
case 3:
//type==3 双余弦子波
parameterDoubleCos DoubleCos_para;
DoubleCos_para.HighCut = 80;
DoubleCos_para.HighPass = 60;
DoubleCos_para.LowCut = 8;
DoubleCos_para.LowPass = 12;
int nK = getK(WaveletLen);
int SpecLen = (int)pow(2.0, nK - 1);
float FreqInc = 1000.0 / (2.0 * WaveletInc * SpecLen);
int ThickRatio = 8;
int DoubleCosSpecLen = int(SpecLen * 2 * ThickRatio);
float *DoubleCosFreq = NULL;
float *DoubleCosSpec = NULL;
DoubleCosFreq = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCosSpec = (float*)calloc(DoubleCosSpecLen, sizeof(float));
DoubleCos_Spec(DoubleCosSpecLen, FreqInc / ThickRatio, DoubleCos_para.LowCut, DoubleCos_para.LowPass,
DoubleCos_para.HighPass, DoubleCos_para.HighCut, DoubleCosFreq, DoubleCosSpec, WaveletAmp);
DoubleCos_Wavelet(DoubleCosSpec, DoubleCosSpecLen, nK, ThickRatio, WaveletLen / 2, WaveletLen, WaveletInc / 1000.0,
WaveletTime, WaveletValue, 1.0);
free(DoubleCosFreq);
free(DoubleCosSpec);
break;
}
float *Signal = NULL;
int lenSignal = WaveletLen + LenRcs - 1;
Signal = (float*)calloc(lenSignal, sizeof(float));
convolve_cwp(LenRcs, 0, Rcs, WaveletLen, 0, WaveletValue, lenSignal, 0, Signal);
memcpy(Syn, Signal + (WaveletLen - 1) / 2, (LenRcs) * sizeof(float));
free(WaveletValue);
free(WaveletTime);
free(Signal);
}
VI main.cpp
#include "ImpToSynProfile.h"
#include"segy.h"
int main(){
const char *filenameInput = "test.segy";
const char *filenameOutput = "Output66.segy";
bhed fileheader;
segy traceheader;
unsigned int nt = 0;
unsigned int sizefileheader = sizeof(fileheader);
unsigned int sizetraceheader = sizeof(traceheader);
unsigned int ncdp = 0;
long long size_file = 0;
long long size_trace = 0;
FILE *fpinput = NULL;
FILE *fpoutput = NULL;
float *dataInput = NULL;
float *dataOutput = NULL;
fpinput = fopen(filenameInput, "rb");
fpoutput = fopen(filenameOutput, "wb");
if (fpinput == NULL) {
printf("can not open %s file\n", filenameInput);
return false;
}
if (fpoutput == NULL) {
printf("can not open %s file\n", filenameOutput);
return false;
}
fread(&fileheader, sizefileheader, 1, fpinput);
nt = fileheader.hns;
_fseeki64(fpinput, 0, SEEK_END);
size_file = _ftelli64(fpinput);
size_trace = nt * sizeof(float) + sizetraceheader;
ncdp = (size_file - (long long)sizefileheader) / size_trace;
_fseeki64(fpinput, sizefileheader, SEEK_SET);
fwrite(&fileheader, sizefileheader, 1, fpoutput);
dataInput = (float*)calloc(nt, sizeof(float));
int WaveletLenTemp = 128;
int WaveletInc = 1;
int WaveletType = 1; //选择子波类型
//遍历每一道
for (int itrace = 0; itrace < ncdp; itrace++)
{
fread(&traceheader, sizetraceheader, 1, fpinput);
fread(dataInput, nt * sizeof(float), 1, fpinput);
//预处理
int *start_end = NULL;
start_end = (int*)calloc(2, sizeof(int));
Preprocessing(dataInput, nt, start_end);
int Start = 0;
int End = 0;
int LenImp = 0;
Start = start_end[0]; //获取有值的部分的起止位置
End = start_end[1];
LenImp = End - Start + 1;
int Boundary = 3; //边界
int LenImp_ex = LenImp + Boundary * 2;
float *Imp = NULL;
float *Rcs = NULL;
float *syn = NULL;
//延展边界2个点
Imp = (float*)calloc(LenImp_ex , sizeof(float));
Rcs = (float*)calloc(LenImp_ex - 1, sizeof(float));
syn = (float*)calloc(LenImp_ex - 1, sizeof(float));
for (int i = Boundary; i < LenImp+Boundary; i++) {
Imp[i] = dataInput[Start-Boundary + i];
}
//上边界
for (int i = 0; i < Boundary; i++) {
Imp[i] = dataInput[Start];
}
//下边界
for (int i = LenImp+Boundary; i < LenImp + Boundary * 2; i++) {
Imp[i] = dataInput[End];
}
//第一步，计算反射系数
Compute_Rcs(Imp,Rcs,LenImp);
//第二步，反射系数褶积子波
Rcs2Syn(Rcs,syn, LenImp_ex - 1, WaveletType, WaveletLenTemp, WaveletInc);
dataOutput = (float*)calloc(nt, sizeof(float));
for (int i = 0; i < LenImp; i++) {
dataOutput[Start + i] = syn[i + Boundary];
// if ((itrace==250)&&(i==250))
// {
// printf("%f", dataOutput[Start+i]);
// }
}
fwrite(&traceheader, sizetraceheader, 1, fpoutput);
fwrite(dataOutput, nt * sizeof(float), 1, fpoutput);
free(start_end);
free(Imp);
free(Rcs);
free(syn);
}
free(dataInput);
free(dataOutput);
fclose(fpinput);
fclose(fpoutput);
return 0;
}