最近写计算程序经常要用到复数运算，以前都是用 gsl 库中对复数运算的支持。这两天看了看C99 标准，发现C99 中对复数运算的支持还是蛮给力的，这里做个总结。C99 中引入了两个关键字 _Complex和 _Imaginary，并且定义了如下的复数类型：

float _Complex
double _Complex
long double _Complex
float _Imaginary
double _Imaginary
long double _Imaginary
这六种复数类型的含义显而易见，这里不多解释。需要特殊指出的是 _Imaginary相关数据类型在当前版本的 gcc 中还没有被支持。因为我现在电脑上没有安装其他的编译器，所以别的编译器能否支持也不清楚。不过 _Imaginary相关数据类型似乎用处不大，至少我没发现什么地方有这个需要。
下面给个简单的程序例子：
#include <stdio.h>

extern double creal(double _Complex);
extern double cimag(double _Complex);
extern double _Complex csqrt(double _Complex);
extern double _Complex csin(double _Complex);

int main(void)
{
    double _Complex a = 1.0 + 2.0 * csqrt(-1);
    double _Complex  b = 5.0 + 4.0 * csqrt(-1);
    printf("a = %f + %fi\n", creal(a), cimag(a));
    printf("b = %f + %fi\n", creal(b), cimag(b));
    a *= b;    printf("a * b = %f + %fi\n", creal(a), cimag(a));
    a = csin(b);    printf("sin(b) = %f + %fi\n", creal(a), cimag(a));
   return 0;
}
计算结果如下：a = 1.000000 + 2.000000i

b = 5.000000 + 4.000000i

a * b = -3.000000 + 14.000000i

sin(b) = -26.186527 + 7.741118i
上面程序中用到了几个函数：
creal() 获得复数的实部

cimag() 获得复数的虚部

csqrt() 计算复数的平方根

csin() 计算复数的sin值
因为程序中没有包含复数支持的头文件，因此需要自己声明这几个函数。我之所以这样做是为了强调 _Complex 是编译器支持的原生类型，不是在头文件中声明的结构体。
C99 标准中之所以使用 _Complex和 _Imaginary而不直接用 complex 和 imaginary 作为关键字，是为了避免与一些现存的C 程序中的变量名发生冲突。如果能够确定自己的代码中没有用到 complex、imaginary 、I 作为变量名，那么可以包含 complex.h 头文件，包含这个头文件后就可以用 complex 代替 _Complex 了。这里需要特别提一句，C99 标准中似乎没有约定单位虚数i 在程序中如何简便的表示。gcc的实现方式是定义了两个宏，_Complex_I和
 I 来表示单位虚数 i 。
因此上面的程序可以简化成如下的样子：
#include <stdio.h>
#include <complex.h>

int main(void)
{
    double _Complex a = 1.0 + 2.0 * I;
    double _Complex  b = 5.0 + 4.0 * I;
    printf("a = %f + %fi\n", creal(a), cimag(a));
    printf("b = %f + %fi\n", creal(b), cimag(b));
    a *= b;    printf("a * b = %f + %fi\n", creal(a), cimag(a));
    a = csin(b);    printf("sin(b) = %f + %fi\n", creal(a), cimag(a));
   return 0;
}


除了这里用到的几个函数之外，gcc 提供的支持复数运算的函数还很多，这里总结如下。

函数名        简介

cacos()        反余弦

casin()        反正弦

catan()        反正切

ccos()        余弦

csin()        正弦

ctan()        正切

cacosh()    反双曲余弦

casinh()    反双曲正弦

catanh()    反双曲正切

ccosh()        双曲余弦

csinh()        双曲正弦

ctanh()        双曲正切

cexp()        指数函数

clog()        对数函数

cpow(x,y)    x 的y 次方

csqrt        平方根

cabs        复数的模

carg        角

conj        共轭

cproj        投影

cimag        取虚部

creal        取实部

参考摘录：https://www.zhihu.com/question/46877027/answer/542742130

前言

高中时候都学过复数，复数表示为a+bi，其中a为实部，b为虚部，i为虚数单位，i^2= -1,然后复数可以用x轴代表实数部，y轴代表虚部的复平面坐标系来表示，最后就是学各种复数运算。

但是什么是虚数？为什么虚数i^2= -1？复数又是怎么来的？什么是复平面？这些都没有解释，所以应用的时候有点小迷茫，最近研究四元数的时候涉及到复数相关的知识，就顺便把以前的这些疑问整理了一下。

 

虚数

虚数就是平方为-1的数，这是我们无法理解的，因为我们以往接触的数学知识里，不管是正数还是负数的平方都是正数，是不存在负数的，这其实是意识形态的问题，就像在早期数学里是没有负数的，数还只停留在可见可数的0、1、2、3....这种，有就有没有就没有，何来负数之说，后来发现有一些实际行为上的计算可以用负数来描述，这才引入了负数。

同样在这里我们不妨假设虚数存在，那么来求解i^2= -1，可以看做是1*i*i= -1，既然按照传统的乘法理解，实数轴上的1（A点）乘了两次一样的数，变成了符号相反的-1（B点）这样难以接受，我们不如将乘法理解成一种转换，实数轴上的1（A点）通过怎样的相同的两次转换能变成-1（B点）呢，答案就是旋转180度（正或负），每次转换就是旋转90度（正或负），我们且以逆时针旋转90度，代表一次转换即乘以虚数i。



 

复数、复数坐标系和复平面

当实数轴上的A点在变换到B点的过程中，首次旋转90度到达C点，坐标为1*i= i， 此时C点已经不在实数轴上，无法用之前的实数轴来描述，此时想要描述这个虚数就必须有一条类似于实数轴的虚数轴来描述它，因为C点的坐标为i，单位纯虚数，也就是说C点和原点所在的这条直线就是虚数的轴，以此便做出了虚数轴。

我们做出的虚数轴与之前存在的实数轴相交于原点，构成一个二维坐标系，此坐标系中落在横轴上的点为实数，落在纵轴上的点为虚数，那么落在坐标轴之外的其他点呢，我们发现，他们既有实数部又有虚数部，这些点就叫做复数，表示为a+bi，这个由实数轴和虚数轴构成的坐标系就叫做复数坐标系，这个复数坐标系所在的平面就叫做复平面。



 

 

目录：

一、简介

二、复数基本性质运算

三、复数二元运算

四、复数三角函数运算

五、C语言复数运算

 

------------------------------------------------------------------------------------------------------

一、简介

在电路设计中我们会遇到一些诸如交流信号过电容、电感，阻抗计算等情况，可以借助Excel中复数计算函数。






 

------------------------------------------------------------------------------------------------------

二、复数基本性质运算

=COMPLEX(-0.5,0.87,"i")创建复数-0.5+0.87i。此处的i在电子中常用作j，可用=COMPLEX(-0.5,0.87,"j")。



 

=IMAGINARY(A2)求出复数的虚部，即虚数i前面的系数。

=IMREAL(A2)求出复数的实部，即实数部分。

 

=IMABS(A2)求出模，即复数的绝对值。注意：模由公式√(x²+y²)得到，其中x实部，y虚部。

=IMARGUMENT(A2)求出复数的辐角，注意是弧度值。

=DEGREES(IMARGUMENT(A2))将弧度转换成角度。






 

=IMCONJUGATE(A2)求出共轭复数。

 

------------------------------------------------------------------------------------------------------

三、复数二元运算

注：对包含公式的单元格引用若报错，参见“Excel高级使用技巧汇总之18、绝对引用与相对引用”

 

=IMSUM(A2,B2)得到复数的和。

注：两个复数虚数字母风格的书写要一致。

例如：(1+2i) +(3-4j) 在Excel中是不规范的，计算不出来结果。



 

=IMSUB(A2,B2)求出复数的差。

=IMPRODUCT(A2,B2)求出复数的乘积。



注：IMPRODUCT支持多个复数之间相乘（最多支持255个）

 

=IMDIV(A2,B2)求出两个复数之间的商。



 

=IMPOWER(A2,1)求出复数的整数幂。

注意，后面的参数可以是复数，即求出复数的倒数

=IMEXP(A2)求出复数的指数（e为底）

=IMLN(A2)求出复数的自然对数。

类似地，可以用对数函数公式

=IMLOG10(A2)

=IMLOG2(A2)

分别求出10为底的对数，2为底的对数。

=IMSQRT(A2)求出复数的平方根。

其实，也可以用Excel公式=IMPOWER(A2,0.5)得到同样的计算结果，换句话说Excel的IMPOWER公式，也支持计算复数的小数形式的幂。

 

------------------------------------------------------------------------------------------------------

四、复数三角函数运算

=IMSIN(A2)求出复数的正弦值。

=IMCOS(A2)求出复数的余弦值。

如果需要正切（正弦函数除以余弦函数）、余切（正切函数的倒数）

另外，像双曲函数、正割、余割函数也可以求解。比较遗憾的是Excel不支持复数求反三角函数值（比如arctan）。

 

------------------------------------------------------------------------------------------------------

五、C语言复数运算

#include "stdio.h"

#include "complex.h"  //加入此头文件需要C99模式支持，详见“Keil编译软件的使用汇总之十六、支持C99模式”

int main(void)

{

   complex double a = 1.0 + 2.0*I, b = 3.0 + 4.0*I; //a = 1 + 2i, b = 3 + 4i

   complex double s_1 = a + b; //加

   complex double s_2 = a - b; //减

   complex double s_3 = a * b; //乘

   complex double s_4 = a / b; //除

 

   printf("s_1的实部为%f, 虚部为%f, 模为%f\n", creal(s_1), cimag(s_1), cabs(s_1));

   printf("s_2的实部为%f, 虚部为%f, 模为%f\n", creal(s_2), cimag(s_2), cabs(s_2));

   printf("s_3的实部为%f, 虚部为%f, 模为%f\n", creal(s_3), cimag(s_3), cabs(s_3));

   printf("s_4的实部为%f, 虚部为%f, 模为%f\n", creal(s_4), cimag(s_4), cabs(s_4));

}

 

complex.h中所包含的常用函数除了基本的加减乘除之外，complex.h中还包含以下函数：

复数的实部creal()

复数的虚部cimag()

模cabs()

相位carg()

指数cexp()

自然对数clog()

幂cpow()

平方根csqrt()

复共轭cong()

正弦csin()

余弦ccos()

正切ctan()

反正弦casin()

反余弦cacos()

反正切catan()

双曲正弦csinh()

双曲余弦ccosh()

双曲正切ctanh()

反双曲正弦casinh()

反双曲余弦cacosh()

反双曲正切catanh()

说明：以上函数默认为双精度版本，在函数名后加f(即float)表示单精度版本，加l(即long)表示长双精度版本。

 

------------------------------------------------------------------------------------------------------

	GDAL的C#版本虽然在很多算法接口没有导出，但是在读写数据中的接口基本上都是完全导出了。使用ReadRaster和WriteRaster方法来进行读写，同时对这两个方法进行了重载，对于常用的数据类型可以不用指定数据类型直接进行读取即可。但是对于复数类型就有点复杂了。下面就针对GDAL如何来读取复数数据来进行一个简单的说明。
	我们知道，在使用GDAL读取数据的时候使用的是ReadRaster这个函数，这个函数重载了6个，函数声明分别如下，以Dataset的ReadRaster为例，Band类中的ReadRaster与之类似。
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Byte[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Int16[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Single[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.Double[], System.Int32, System.Int32, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
OSGeo.GDAL.Dataset.ReadRaster(System.Int32, System.Int32, System.Int32, System.Int32, System.IntPtr, System.Int32, System.Int32, OSGeo.GDAL.DataType, System.Int32, System.Int32[], System.Int32, System.Int32, System.Int32);
	前五个函数的声明基本上完全一致，除了第五个参数，分别是Byte[]，Int16[]，Int32[]，Single[]，Double[]。分别对应的是将图像中的像素值使用8U，16U\16S，32U\32S，32F和64F。如果图像数据是普通的数据（这里指的是除复数图像之外），那么这五个函数完全可以满足所有的需求。但是如果图像是复数图像（复数图像就是图像的像素值是由复数组成的，比如通过傅立叶变换后的图像），这时，上面的五个函数就有点爱莫能助了，就需要上面的第六个函数来出场救急了。
	第六个函数的声明和C\C++的声明完全一致，使用C\C++的同学对这个接口应该比较熟悉。这个函数可以替代上面五个函数来使用，但是比较麻烦的时，就是传入的参数会比较多，首先看一个例子。使用普通的方式来读取一个单波段8U的数据（Landsat5的一个数据）。从数据中读取100×100大小的数据，直接使用ReadRaster的Byte接口。
static void ReadRaster()
{
	string strFile = @"D:\Data\landsat\etp139r26_5t19900902\p139r26_5t19900902_nn1.tif";

	OSGeo.GDAL.Gdal.AllRegister();
	OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(strFile, OSGeo.GDAL.Access.GA_ReadOnly);
   
	int []bandmap = new int [1];
	bandmap[0] = 1;

	//读取Byte类型
	Byte[] data = new Byte[100 * 100];
	OSGeo.GDAL.CPLErr err = ds.ReadRaster(3000, 2000, 100, 100, data, 100, 100, 1, bandmap, 0, 0, 0);
	
	//将读取的数据输出
	for (int i = 0; i < 100; i++ )
	{
		for (int j = 0; j < 100; j++)
		{
			Console.WriteLine("{0}", ddata[i*100+j]);
		}
	}
	Console.WriteLine("done\n");
}
	下面看看如何使用第六个接口来完成上面的功能，代码如下：
static void ReadRaster()
{
	string strFile = @"D:\Data\landsat\etp139r26_5t19900902\p139r26_5t19900902_nn1.tif";

	OSGeo.GDAL.Gdal.AllRegister();
	OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(strFile, OSGeo.GDAL.Access.GA_ReadOnly);
   
	int []bandmap = new int [1];
	bandmap[0] = 1;

	//读取Byte类型
	Byte[] data = new Byte[100 * 100];
	
	//构造一个IntPtr，大小为100×100个Byte
	IntPtr ptr = Marshal.AllocCoTaskMem(data.Length);
	//需要指明读取的数据类型
	OSGeo.GDAL.CPLErr err = ds.ReadRaster(3000, 2000, 100, 100, ptr, 100, 100, OSGeo.GDAL.DataType.GDT_Byte, 1, bandmap, 0, 0, 0);
	//将读取到的数据拷贝到data中
	Marshal.Copy(ptr, data, 0, data.Length);

	//将读取的数据输出
	for (int i = 0; i < 100; i++ )
	{
		for (int j = 0; j < 100; j++)
		{
			Console.WriteLine("{0}", ddata[i*100+j]);
		}
	}
	Console.WriteLine("done\n");
}
	从上面的代码中可以看出，使用第六个函数的时候，需要一个IntPtr的类型，使用Marshal类中的AllocCoTaskMem函数来进行分配内存空间（不知道C#中的说法是啥，按照C\C++中应该就是分配空间）。需要注意的是，分配的大小是按照字节（Byte）为单位的。接下来在调用ReadRaster时需要指定数据类型，调用之后，像素值应该都存储在这个IntPtr中了，然后再使用Marshal类中的Copy函数将IntPtr中的数据拷贝到Byte数组data中即可。
	通过上面的代码，我们大致知道第六个ReadRaster的大致用法，接下来我们通过这个方法来读取一个复数类型的图像。没有复数类型的图像，可以使用ENVI随便打开一个图像，然后在Transform菜单下有个FFT，将打开的数据进行傅立叶变换，输出的结果就是一个复数图像。这里我将上面使用的Landsat的数据使用FFT转了下。使用下面的代码进行打开。
static void ReadRaster()
{
	string strFile = @"D:\Data\landsat\etp139r26_5t19900902\p139r26_5t19900902_nn1_fft.tif";

	OSGeo.GDAL.Gdal.AllRegister();
	OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(strFile, OSGeo.GDAL.Access.GA_ReadOnly);
   
	int []bandmap = new int [1];
	bandmap[0] = 1;


	//读取复数类型，由于一个复数由两个数据组成，所以读取100×100的像元需要2倍普通图像的大小
	double[] data = new double[100 * 100 * 2];
	
	//构造一个IntPtr，大小为100×100×2个double，一个double为4个byte，所以下面的长度要乘以8
	IntPtr ptr = Marshal.AllocCoTaskMem(data.Length*8);
	//需要指明读取的数据类型
	OSGeo.GDAL.CPLErr err = ds.ReadRaster(0, 0, 100, 100, ptr, 100, 100, OSGeo.GDAL.DataType.GDT_CFloat64, 1, bandmap, 0, 0, 0);
	//将读取到的数据拷贝到data中
	Marshal.Copy(ptr, data, 0, data.Length*8);

	//将读取的数据输出
	int in1 = 0;
	for (int i = 0; i < 100; i++ )
	{
		for (int j = 0; j < 100; j++)
		{
			Console.WriteLine("({0}*i+{1})", ddata[in1++], ddata[in1++]);
		}
	}
	Console.WriteLine("done\n");
}
	在上面的代码中，一个复数由实部和虚部组成，所以一个像素值就对应于普通图像的2个像素值，故读取100×100大小的数据就需要分配普通图像2倍的大小。除此之外与读取普通图像一样，不过在最后获取像素值的时候，复数图像的像素值是顺序存储的，即实部,虚部,实部,虚部...这样的顺序来进行存储。上面的代码运行的前五个像素值如下图所示：

	通过ENVI软件查看该图像的前五个像素值如下图所示（图中红色框中就是像素值），对比出来可以看到读取出来的实部和虚部完全一样。