C/C++实现合成地震记录

来源：互联网收集：自由互联发布时间：2022-06-23

C/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;
}

上一篇：Flink状态管理
下一篇：没有了

网友评论

相关栏目

C/C++实现合成地震记录

相关文章