CSAMT可控源音频电磁法仪器研制与离散富里叶变换(DFT).doc

CSAMT可控源音频电磁法仪器研制与离散富里叶变换(DFT)*离散富里叶变换(DFT)在地质仪器中应用 对于大地系统的分析，电法勘探中主要是用电流作为输入信号，电压作为输出信号，大地是信号传递系统，如下图所示： Vt发送机输出电流时域信号方波I (t)频域信号I (f)接收机获得信号时域，准方波大地系统时域脉冲响应函数 h (t)频域传递函数H (f) I If Vf fo f fo f 图（1） 上图中的It是时域电流的波形，If是频域中电流的频谱。Vt是时域的电压波形，Vf是频域中电压的频谱。 以上时域信号都是可以测量获得的，测量的办法时对电流和电压的波形进行采样和模数转换（A/D），如下图所示： (I) t h (K0) h (K1023) (V) t K0 K1023 t T N 图（2） 在电法勘探中使用的波形是连续的周期性的稳态波形。这是为了用频率域分析信号的要求，要求周期性是为了用计算机使用的离散采样信号来分析物理的连续信号，用傅立叶级数计算来代替无穷时间的连续积分。 对于时间域的信号It输入到大地的系统，经过大地系统的作用输出的时间域信号Vt。Vt与It之间的关系是： Vt = It（卷积）h(t) (1)这是一个积分方程，很难用计算机计算处理。 而在频率域中Vf和If的关系是一个简单的代数式： Vf = If * H（f） （2） 我们的目的是要求得大地传递函数H（f），即大地的电阻率频谱函数，这里的H（f）传递函数是一个复数，它可以写成用振幅 | H（f）| 和相位R 表示的e指数形式： j*R H (f) = | H（f）| * e （3） 上式中j是虚部-1 利用公式（2）和（3）可得到下式： j*V Vf |Vf | * e |Vf | j*(V -I ) H (f) = = = * e (4) j*IIf |If | * e |If |上式中 |Vf | 是输入电压的频谱幅值，V是输入电压的相位 |If |是输入电流的频谱幅值，I是输入电压的相位利用公式（3）和公式（4）可得到下式：j*R |Vf | j*(V -I ) | H（f）| * e = * e （5） |If |根据（5）式我们可得到大地系统的复电阻幅值： |Vf | | H（f）| = （6） |If |大地系统的相位： R =V I （7）到此我们意义竟清楚地看出，在频率域中求大地系统的复电阻的幅值和相位用到的计算方式都是简单的代数方法，计算机很容易计算处理。但问题的关键是我们如何求出|Vf |，V，|If |，I这些量。这正是我们要在下文中讨论的DFT离散傅立叶变换算法。用高速A/D对电流和电压的时域信号进行一个完整周期的采样，然后分别对两个数据序列进行离散傅立叶变换，就能用加法和乘法来算出|Vf |，V，|If |，I 为了说明离散的傅立叶变换，我们先要从速读的傅立叶变换说起，傅立叶积分定义为： -j2ftH（f）= h(t) * e dt (8) -上式中H（f）是频率域内的频谱函数 h(t)是时间域内的时间波形信号函数因为 -j2fte = cos(2ft)- jsin(2ft) (9)用（8）式和（9）式可得到以三角函数表达的傅立叶变换式： H（f）= h(t) * cos(2ft)- jsin(2ft) dt （10） - 对于周期性的时域信号h(t)可以用下面傅立叶级数形式来代替上式： N-1 H（f n）= h(k) * cos(2f n * k)- jsin(2f n * k) K=0 N-1 N-1= h(k) * cos(2f n * k)- j * sin(2f n * k) K=0 K=0= Af n j* B f n (11) 上式中Af n是频谱函数H（f n）的f n谐波分量的实数部分系数 B f n是频谱函数H（f n）的f n谐波分量的虚数部分系数 如果f n是基波，则f 0=1，上式变得更加简单。我们在地质仪器中为了提高精度只对方波的基波进行傅立叶变换，实践证明对基波采1024个样，A/D转换分辨达到12bit时，可以得到幅度和相位分辨率0.1%，相位精度在2个毫弧度左右。 我们得到对方波信号的基波离散傅立叶变换的公式如下：1023 2 2H（f0）= h(k) * cos( * k)- j * h(k)sin( * k) （12） K=0 1024 1024 参看图（2）上式中，h(k)是波形中第K个采样点的A/D转换得的数据，它是一个二字节 2的整数，让 =，上式变为： 10241023 H（f0）= h(k) * cos( * k)- j * h(k)sin( * k) K=0 =A f0 j* B f0 （13）上式中的 1023 A f0 = h(k) * cos( * k) （14） K=0 1023 B f0 = h(k) * sin( * k) （15） K=0 A f0是频谱函数的实部，B f0是频谱函数的虚部 用e指数形式，（13）式可写成 B f0 jarctg 2 2 A f0H（f0）= （A f0） + （B f0） * e j*f0 = | H（f0）| * e (16) 2 2 上式中 |H（f0）|= （A f0） + （B f0） （17） 它是基波f0的幅值 B f0 f0 = arctg A f0 这是基波f0 的相位。注意这个相位是随机的，但是对电流I波形和电压V波形进行同步采样，用傅立叶变换求出的I电流相位）和V（电压相位），虽然它们都是随机的，随机的值取决与开始采样的那时刻，但V I = 电压对电流的相移是固定的，它是由大地系统的传递函数决定的。它就是我们要测量的物理参数，它就是复电阻率的相位。综上所述，要用傅立叶变换对大地系统测量频谱特性的方法步骤如下： 1：让发送机想地下供电，发送机输出的电流波形是方波，频率是精确的（1/10万 精确度），电流的幅度不要求很精确（5%精确度），用接收机中A/D转换器对方波采样1024个数据，A/D的分辨率为12bit。为了使放大器部分与发送机部分的输出端上的高压隔离，需要一个隔离电压在2500V以上，频率响应在100KHZ以上的隔离放大器。发送机的电流值用0.1的精密电阻（自制的锰铜丝电阻）转换为电压值送给隔离变压器，隔离放大器输出的电压与发送机输出端高压隔离，剩下的只是I*0.1的电压VI值。这个电压就是发送机输出的电流波形，如下图所示：高压隔离变压器放大器A/D发送机 A微处理器 B 接收机 2:接收机对发送机输出电流波形采样的同时,同步地对来自大地某测点的V波形进行采样1024个点.A/D的分辨率是12bit,为了提高测量精度对来自大地M,N接收电极上的信号要进行放大,放大到A/D能接收的最高电压的一半以上.放大器的放大倍数是可以用计算机程序控制的(1100),这样就形成了一个浮点A/D转换器,它的分辨率是12bit (1/4096),而电压测量的范围是满度满度值: + 50mV + 5V最高分辨率: 50 mV /2048 = 24.4 25V 对于野外信号有25V的分辨率是足够了,再高会被电极的极差漂移,电线上的干扰,放大器的噪声淹没. 3:分别用离散傅立叶变换酸楚电流波形的频谱幅值和相位得到: jII(f0)= |I(f0)| * e jVV(f0)= |V(f0)| * e 4:把电压的幅值除以电流的幅值得到地质体的复电阻幅值. |V(f0)| R(f0)= |I(f0)| 5:把电压的相位减去电流的相位,得到地质体的复电阻的相位 R = V I 6:把复电阻幅值乘上野外布极的装置系数K值,得到地质的视复电阻率s = K * | R(f0)| 7:改变发送机的输出电流的频率后再重复以上6个步骤,这样就可以得到地质体对不同频率的响应,得到完整的地质体复电阻率的幅值随频率而变化的相位图,如下图所示: H, 幅值曲线 相位曲线 f 该图可在野外接收机上显示出来 8:以上步骤还缺少一个仪器自身放大器的频谱特性对测量带来的影响问题,本文用另一篇文章来说明校正方法./*/* 离散富里变换程序 */*/*DFT()函数中计算sin(M*2PI/1024)使用查表法*/*存放在CORD(ROM)中的sin()表*/static float code sin_table_float260= 0.0000000, 0.0061359, 0.0122715, 0.0184067, 0.0245412, 0.0306748, 0.0368072, 0.0429383, 0.0490677, 0.0551952, 0.0613207, 0.0674439, 0.0735646, 0.0796824, 0.0857973, 0.0919090, 0.0980171, 0.1041216, 0.1102222, 0.1163186, 0.1224107, 0.1284981, 0.1345807, 0.1406582, 0.1467305, 0.1527972, 0.1588582, 0.1649131, 0.1709619, 0.1770042, 0.1830399, 0.1890687, 0.1950903, 0.2011046, 0.2071114, 0.2131103, 0.2191012, 0.2250839, 0.2310581, 0.2370236, 0.2429802, 0.2489276, 0.2548656, 0.2607941, 0.2667128, 0.2726214, 0.2785197, 0.2844075, 0.2902847, 0.2961509, 0.3020059, 0.3078496, 0.3136818, 0.3195020, 0.3253103, 0.3311063, 0.3368899, 0.3426607, 0.3484187, 0.3541635, 0.3598951, 0.3656130, 0.3713172, 0.3770074, 0.3826834, 0.3883450, 0.3939920, 0.3996242, 0.4052413, 0.4108432, 0.4164295, 0.4220003, 0.4275551, 0.4330938, 0.4386162, 0.4441221, 0.4496113, 0.4550836, 0.4605387, 0.4659765, 0.4713967, 0.4767992, 0.4821838, 0.4875502, 0.4928982, 0.4982277, 0.5035384, 0.5088301, 0.5141028, 0.5193560, 0.5245897, 0.5298036, 0.5349976, 0.5401714, 0.5453250, 0.5504580, 0.5555702, 0.5606616, 0.5657318, 0.5707808, 0.5758082, 0.5808139, 0.5857978, 0.5907597, 0.5956993, 0.6006165, 0.6055111, 0.6103828, 0.6152316, 0.6200572, 0.6248595, 0.6296383, 0.6343933, 0.6391245, 0.6438316, 0.6485144, 0.6531729, 0.6578067, 0.6624158, 0.6669999, 0.6715590, 0.6760927, 0.6806010, 0.6850837, 0.6895406, 0.6939715, 0.6983762, 0.7027547, 0.7071068, 0.7114322, 0.7157308, 0.7200025, 0.7242471, 0.7284644, 0.7326543, 0.7368166, 0.7409511, 0.7450578, 0.7491364, 0.7531868, 0.7572088, 0.7612024, 0.7651672, 0.7691033, 0.7730104, 0.7768885, 0.7807372, 0.7845566, 0.7883464, 0.7921066, 0.7958369, 0.7995372, 0.8032075, 0.8068476, 0.8104572, 0.8140363, 0.8175848, 0.8211025, 0.8245893, 0.8280451, 0.8314696, 0.8348629, 0.8382247, 0.8415550, 0.8448536, 0.8481203, 0.8513552, 0.8545580, 0.8577286, 0.8608670, 0.8639728, 0.8670462, 0.8700870, 0.8730950, 0.8760701, 0.8790122, 0.8819213, 0.8847971, 0.8876396, 0.8904487, 0.8932243, 0.8959662, 0.8986745, 0.9013488, 0.9039893, 0.9065957, 0.9091680, 0.9117060, 0.9142098, 0.9166791, 0.9191139, 0.9215140, 0.9238795, 0.9262102, 0.9285061, 0.9307669, 0.9329928, 0.9351835, 0.9373390, 0.9394592, 0.9415441, 0.9435934, 0.9456073, 0.9475856, 0.9495282, 0.9514350, 0.9533060, 0.9551412, 0.9569404, 0.9587035, 0.9604305, 0.9621214, 0.9637761, 0.9653944, 0.9669765, 0.9685221, 0.9700313, 0.9715039, 0.9729400, 0.9743394, 0.9757021, 0.9770281, 0.9783174, 0.9795698, 0.9807853, 0.9819639, 0.9831055, 0.9842101, 0.9852777, 0.9863081, 0.9873014, 0.9882576, 0.9891765, 0.9900582, 0.9909027, 0.9917098, 0.9924796, 0.9932119, 0.9939070, 0.9945646, 0.9951847, 0.9957674, 0.9963126, 0.9968203, 0.9972904, 0.9977230, 0.9981181, 0.9984756, 0.9987954, 0.9990777, 0.9993224, 0.9995294, 0.9996988, 0.9998306, 0.9999247, 0.9999812, 1.0000000, 0.9999812, 0.9999247, 0.9998306, ;/* cosine table */static float code cos_table_float260= 1.0000000, 0.9999812, 0.9999247, 0.9998306, 0.9996988, 0.9995294, 0.9993224, 0.9990777, 0.9987954, 0.9984756, 0.9981181, 0.9977230, 0.9972904, 0.9968203, 0.9963126, 0.9957674, 0.9951847, 0.9945646, 0.9939070, 0.9932119, 0.9924796, 0.9917098, 0.9909027, 0.9900582, 0.9891765, 0.9882576, 0.9873014, 0.9863081, 0.9852777, 0.9842101, 0.9831055, 0.9819639, 0.9807853, 0.9795698, 0.9783174, 0.9770281, 0.9757021, 0.9743394, 0.9729400, 0.9715039, 0.9700313, 0.9685221, 0.9669765, 0.9653944, 0.9637761, 0.9621214, 0.9604305, 0.9587035, 0.9569404, 0.9551412, 0.9533060, 0.9514350, 0.9495282, 0.9475856, 0.9456073, 0.9435934, 0.9415441, 0.9394592, 0.9373390, 0.9351835, 0.9329928, 0.9307669, 0.9285061, 0.9262102, 0.9238795, 0.9215140, 0.9191139, 0.9166791, 0.9142098, 0.9117060, 0.9091680, 0.9065957, 0.9039893, 0.9013488, 0.8986745, 0.8959662, 0.8932243, 0.8904487, 0.8876396, 0.8847971, 0.8819213, 0.8790122, 0.8760701, 0.8730950, 0.8700870, 0.8670462, 0.8639728, 0.8608670, 0.8577286, 0.8545580, 0.8513552, 0.8481203, 0.8448536, 0.8415550, 0.8382247, 0.8348629, 0.8314696, 0.8280451, 0.8245893, 0.8211025, 0.8175848, 0.8140363, 0.8104572, 0.8068476, 0.8032075, 0.7995372, 0.7958369, 0.7921066, 0.7883464, 0.7845566, 0.7807372, 0.7768885, 0.7730104, 0.7691033, 0.7651672, 0.7612024, 0.7572088, 0.7531868, 0.7491364, 0.7450578, 0.7409511, 0.7368166, 0.7326543, 0.7284644, 0.7242471, 0.7200025, 0.7157308, 0.7114322, 0.7071068, 0.7027547, 0.6983762, 0.6939715, 0.6895406, 0.6850837, 0.6806010, 0.6760927, 0.6715590, 0.6669999, 0.6624158, 0.6578067, 0.6531729, 0.6485144, 0.6438316, 0.6391245, 0.6343933, 0.6296383, 0.6248595, 0.6200572, 0.6152316, 0.6103828, 0.6055111, 0.6006165, 0.5956993, 0.5907597, 0.5857978, 0.5808139, 0.5758082, 0.5707808, 0.5657318, 0.5606616, 0.5555702, 0.5504580, 0.5453250, 0.5401714, 0.5349976, 0.5298036, 0.5245897, 0.5193560, 0.5141028, 0.5088301, 0.5035384, 0.4982277, 0.4928982, 0.4875502, 0.4821838, 0.4767992, 0.4713967, 0.4659765, 0.4605387, 0.4550836, 0.4496113, 0.4441221, 0.4386162, 0.4330938, 0.4275551, 0.4220003, 0.4164295, 0.4108432, 0.4052413, 0.3996242, 0.3939920, 0.3883450, 0.3826834, 0.3770074, 0.3713172, 0.3656130, 0.3598951, 0.3541635, 0.3484187, 0.3426607, 0.3368899, 0.3311063, 0.3253103, 0.3195020, 0.3136818, 0.3078496, 0.3020059, 0.2961509, 0.2902847, 0.2844075, 0.2785197, 0.2726214, 0.2667128, 0.2607941, 0.2548656, 0.2489276, 0.2429802, 0.2370236, 0.2310581, 0.2250839, 0.2191012, 0.2131103, 0.2071114, 0.2011046, 0.1950903, 0.1890687, 0.1830399, 0.1770042, 0.1709619, 0.1649131, 0.1588582, 0.1527972, 0.1467305, 0.1406582, 0.1345807, 0.1284981, 0.1224107, 0.1163186, 0.1102222, 0.1041216, 0.0980171, 0.0919090, 0.0857973, 0.0796824, 0.0735646, 0.0674439, 0.0613207, 0.0551952, 0.0490677, 0.0429383, 0.0368072, 0.0306748, 0.0245412, 0.0184067, 0.0122715, 0.0061359, 0.0000000, -0.0061359, -0.0122715, -0.0184067, ;static float Delta_Rad=6.283/1024;uchar ad2048;/*以字节形式存放的1025个a/d采样值,低字在前,高字节右后*/float AD_Samp1024; /*采样所获得的2048个字节,1024个浮点的XDATA数组*/int ADBUF1024;/*当前A/D转换器获得的1024个整数值*/float TOP_FV;/*波形的峰值(浮点数)*/float AVER_FV;/*波形的峰值(浮点数)*/float DLTV1M;/*一次场的面积*/float DLTV2M;/*二次场的面积*/float DLTVM;/*二次场的面积除以一次场的面积得充电率*/long TOP_V;/*波形的峰值(长整数)*/long NTOP_V;/*波形的谷值(长整数)*/long AVER_V;/*波形的平均值(长整数)*/*main()*/main()Set_SP_GAN();/*自电补偿和自动增益设置*/printf(Start Samp 1024n);Samp1024(0.25);/*对频率为0.25hz的波形采12位的a/d值1024个,存放在uchar ad2048中*/ clear_screen();Disp_ad_wav();/*把ad2048中的1024个a/d值显示*/_getkey(); clear_screen();printf(Get_Float_Adn); Get_Float_Ad();/*把ad2048中的偏移码a/d采样值转化成AD_Samp1024的浮点数*/printf(AD_Samp1024 nDFT start !n); DFT();fz=fz/GAN;printf(DFT resulte:nfz=%fnxw=%fn,fz,xw);printf(Pres Any Key Ret);_getkey();/*/* 离散富里叶变换 */* 1024 1024 */* (nf0)=h(m)*cos(2*pi*n*m)/1024)-jh(m)*sin(2*pi*n*m)/1024)/1024 */* m=0 m=0 */*/ float fz,xw;/*fz是DFT计算的复数的幅值,xw是相位弧度值*/void DFT()/* 测试离散富里