声波有限差分法正演模拟c语言程序
C语言中的音频处理方法
C语言中的音频处理方法音频处理在计算机科学领域中扮演着重要的角色，无论是进行语音识别、音乐合成还是语音通信，都需要对音频进行处理。
C语言作为一种高效、灵活的编程语言，提供了许多方法来实现音频处理。
本文将介绍一些常用的C语言音频处理方法。
一、音频采样与表示在音频处理中，我们通常将声波信号转化为数字形式进行处理。
这涉及到音频采样与表示。
C语言中，我们可以使用short或float类型的数组来表示音频数据。
short类型适用于低位数、高速的音频处理，而float类型适用于高位数、低速的音频处理。
根据实际需求选择合适的数据类型。
二、音频录制与播放音频录制与播放是音频处理的基础。
C语言提供了一些库来实现音频录制与播放的功能。
例如，我们可以使用PortAudio库来进行跨平台的音频录制与播放。
使用PortAudio库，我们可以直接处理音频输入和输出流，实现实时音频处理的功能。
另外，还可以使用OpenAL库来实现3D音频的播放效果，通过控制音频源的位置、方向和距离，让听者有身临其境的感觉。
三、音频滤波音频滤波是音频处理中常用的技术，用于消除或增强特定频率的信号。
C语言提供了一些滤波器的实现方法。
例如，我们可以使用FIR（有限脉冲响应）滤波器来实现音频滤波。
FIR滤波器是一种线性相位滤波器，可以通过计算每个样本的加权平均值来实现滤波效果。
此外，我们还可以使用IIR（无限脉冲响应）滤波器，它可以实现更复杂的滤波器特性。
四、音频压缩与解压缩音频压缩与解压缩是音频处理中常用的技术，用于减小音频文件的大小，提高传输效率。
C语言中，我们可以使用库来实现音频压缩与解压缩。
例如，我们可以使用LAME库来实现MP3音频的压缩与解压缩。
LAME库是一个开源的音频编码库，可以提供高质量的音频压缩效果。
此外，还可以使用其他压缩算法，如AAC、OGG等。
五、音频特征提取音频特征提取是音频处理中的重要任务，用于从音频信号中提取出有用的特征信息。
一个简单的dspC语言例子
一个简单的dsp C语言例子开发平台: CCS集成开发环境通过这个简单的例子, 可以大致了解用C语言开发dsp程序的原理。
程序要求: 用C语言编写产生正弦调幅波信号的源程序；正弦调幅波的公式在离散域中的表示:y(n) = (1 + M*sin(2 * PI * fb / fs * n)) * sin(2 * PI * fa / fs * n);编写文件1.sin_am.c#include<stdio.h>#include<math.h>#define TRUE 1#define pi 3.1415926536int y[500],i;float M;void main(){puts("amplitude modulation sinewave example started.\n");M = 50;for(i = 0; i < 500; i++)y[i]= 0;while(TRUE){for(i = 0; i < 500; i++)y[i]=(int)((1 + M / 100 * sin(i * 2 * pi * 20 / 4000))* sin(i * 2 * pi * 200 / 4000)* 16384);puts("program end");}}2.sin_am_v.asm (reset vector file).title "sin_am_v.asm".sect ".vectors".ref _c_int00RESET:B _c_int00.end..3.sin_am.cmdsin_am.objsin_am_v.obj-m sin_am.map-o sin_am.outMEMORY{PAGE 0:EPROG: origin = 0x1400, len = 0x7c00 VECT: origin = 0xff80, len = 0x80PAGE 1:USERREGS: origin = 0x60, len = 0x1c IDATA: origin = 0x80, len = 0x3000 }SECTIONS{.vectors:>VECT PAGE 0.text:>EPROG PAGE 0.cinit:>EPROG PAGE 0.bss:>IDATA PAGE 1.const:>IDATA PAGE 1.switch:>IDATA PAGE 1.system:>IDATA PAGE 1.stack:>IDATA PAGE 1}"*.cmd"文件说明:链接命令文件是实现对段的存储空间位置的定位, C语言程序中常用已初始化和未初始化段如下:已初始化段包括:.init 存放C程序中的变量的初值和常量, 放在ROM和RAM 中均可, 一般属于PAGE 0.const 存放C程序中的字符常量、浮点常量和用const声明的常量, 放在ROM和RAM中均可, 一般属于PAGE 1.text 存放C程序代码, 放在ROM和RAM中均可, 一般属于PAGE 0.switch 存放C程序中的语句的跳针表, 放在ROM和RAM中均可, 一般属于PAGE 0未初始化段包括:.bss 为C程序中的全局和静态变量保留存储空间, 一般存放于RAM中, 属于PAGE 1.stack 为C程序系统堆栈保留存储空间, 用于保存返回地址、函数间的参数传递、存储局部变量和保存中间结果, 一般存放于RAM中, 属于PAGE 1.sysmem 用于C程序中malloc、calloc和realloc函数动态分配存储空间, 一般存放于RAM中, 属于PAGE 14.vary_M.gelmenuitem "Myfunctions"slider vary_M(0, 100, 10, 1, Amount_of_modulation){M = Amount_of_modulation;}该文件用于调试的时候可随意改变变量M的值, 该文件通过file->load GEL File添加到工程中, 调试的时候可选择GEL->My Functions->vary_M来打开vary_M滑动条组件。
C语言音频处理与信号处理
C语言音频处理与信号处理音频处理和信号处理是计算机科学领域中重要的研究方向。
C语言作为一种高级编程语言，在音频处理和信号处理的应用中扮演着重要的角色。
本文将会介绍C语言在音频处理和信号处理方面的应用，并探讨一些常见的技术和算法。
1. 概述1.1 音频处理音频处理是指对音频信号进行各种操作和转换的过程。
它包括音频信号的采集、存储、压缩、解压缩、滤波、降噪、音频合成等。
在音频处理中，C语言提供了丰富的库函数和工具，方便开发者进行各种操作。
1.2 信号处理信号处理是指对信号进行数字化处理的过程。
它包括信号的采集、采样、量化、编码、滤波、频谱分析、时频分析等。
C语言具有较高的效率和灵活性，适用于实时信号处理和非实时信号处理。
2. 音频处理算法2.1 声音采集与播放在C语言中，可以通过使用库函数将声音信号从外部设备（例如麦克风）中采集到计算机，并实现声音的播放。
常用的库函数有"waveInOpen"、"waveInStart"、"waveOutOpen"和"waveOutWrite"等。
2.2 声音压缩与解压缩声音压缩可以将原始的音频信号转化为更小的尺寸，以便节省存储空间和传输带宽。
在C语言中，可以使用各种压缩算法，如MP3、AAC等，通过调用相关的库函数实现。
2.3 声音滤波声音滤波是对音频信号进行频率滤波的过程，以达到去除噪音或者增强特定频率的目的。
在C语言中，可以通过设计数字滤波器并使用库函数来对音频信号进行滤波操作。
2.4 声音合成声音合成是通过对已有的音频片段进行组合和调整，生成新的声音信号。
C语言提供了合成音频的相关库函数，如"mixSounds"和"combineSounds"等，开发者可以利用这些库函数实现声音的合成。
3. 信号处理算法3.1 信号采集与处理C语言提供了许多函数和工具来采集和处理信号。
一种利用有限差分来正演模拟声波波形的方法
一种利用有限差分来正演模拟声波波形的方法
王晓飞;刘海涅
【期刊名称】《科技风》
【年(卷),期】2012(000)022
【摘要】声波波形对于研究并旁地层的情况有着非常重要的意义.有限差分是常用的正演方法.本文利用有限差分的方法对软地层和硬地层的不同模型进行计算,并对结果进行了分析.
【总页数】3页(P42-44)
【作者】王晓飞;刘海涅
【作者单位】中海油田服务股份有限公司油田技术事业部,北京市101149;中海油田服务股份有限公司油田技术事业部,北京市101149
【正文语种】中文
【相关文献】
1.利用哈特莱变换进行井间声波波场正演模拟 [J], 刘迎曦;张霖斌
2.流固边界耦合介质高阶有限差分地震正演模拟方法 [J], 吴国忱;李青阳;吴建鲁;梁展源
3.一种新型有限差分网格剖分方法在大地电磁一维正演中的应用 [J], 张辉;唐新功
4.利用远震波形反演和宽频带地震波正演模拟推断2008年汶川地震的破裂过程[J], Takeshi Nakamur;Seiji Tsuboi;Yoshiyuki Kaneda;Yoshiko Yamanaka;付萍杰;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
5.基于三维有限差分方法的三分量\r感应测井正演模拟 [J], 郭晨;陈晓亮;卢圣鹏
因版权原因，仅展示原文概要，查看原文内容请购买。
声波方程有限差分正演
题目：使用Ricker 子波，刚性边界条件，并且初值为零，在均匀各向同性介质条件下，利用交错网格法求解一阶二维声波方程数值解。
解：一阶二维声波方程：22222221zPx P t P c ∂∂+∂∂=∂∂ (1)将其分解为：21P c t Px P z x z x z V V x z V tV t ∂∂∂⎧=+⎪∂∂∂⎪∂∂⎪=⎨∂∂⎪∂∂⎪=⎪∂∂⎩(2)对分解后的声波方程进行离散，可得到：112211,-1,,,122[]N n n n n m i m j i m j xi j xi j m t VVc P P h +-+---=∆=+-∑ 112211,1,,,122[]Nn n n n m i j m i j m zi j zi j m t VV c P P h +-++---=∆=+-∑ 1111212222,,m 1,,,,11[]Nn n n n n n i ji jmxi j xi m j zi j m zi j m m tc PP cVVVVh+++++++-+--=∆=+-+-∑h z x =∆=∆针对公式(1)，使用二阶中心差商公式：2P(,,1)2(,,)(,,1)i j n P i j n P i j n t +-+-∆222(1,,)2(,,)(1,,)(,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆=⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(3)变形：P(,,1)=2(,,)(,,1)i j n P i j n P i j n +--2222(1,,)2(,,)(1,,)t (,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆+∆⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(4)对离散格式作时间和空间三重Fourier 变换：0P(,,)(,,)x z i j n P k k w ↔ ，0P(,,1)(,,)*exp()x z i j n P k k w iw t +↔∆0P(1,,)(,,)*exp(k )x z x i j n P k k w i x +↔-∆，0z P(,1,)(,,)*exp(k )x z i j n P k k w i z +↔-∆对公式(4)进行Fourier 变换：2222exp()2exp()h exp()2()exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆=--∆+∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦2222exp()2exp()h exp()2()=exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆-+-∆∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦222222sin sin 22sin (2x z k x k zw tt c h∆∆+∆=∆） (5) 公式(5)右端必须满足下列条件：22222sin sin 220(x z k x k zt c h∆∆+≤∆≤）1 取x k 和z k 最大值，即=x x k π∆，z =k z π∆，则有：22220t c h≤∆≤1因此tc ∆≤即为所求得的稳定性条件。
有限差分法地震波传播数值模拟
=
kΔx
=
2π λ
Δx ≤ 1
，即只要一个波长包含几个空间步
长，随着差分精度2M的提高，上述高阶差分解法产生
的数值频散会逐渐减小。
不同差分精度空间频散曲线
不同差分精度时间频散曲线
五、边界问题
自由边界条件
内部边界条件 吸收边界条件
设计吸收边界条件的目标：
z 方程+边界条件数学上是非病态的
连续
问题 z 方程+边界条件可以近似描述无限介质中的物理过程 z 边界条件和内部点的计算方式是相容、不冲突的
-----------J.M. Carcione
地震波传播数值模拟应用领域
地震波传播理论
数据采集
理论指导 物性参数
研究传播规律
正演模拟
指导设计 观测系统
验证
地震解释
提供理论数据 试验处理流程
数据处理
提供正演方法
岩石物理
参数反演
断层下覆界面反射能量强
炮点
T=2000ms
炮点 T=2300ms
炮点位于11km处的单炮记录
?21?4?1?6?1?m?2m?1222m246333m24lllol622m32m?m4?m6??m?m?2m?m?2?c1m??m??c2?m??c3???m??cm??m??1??0????0???m????0??35?1?133335?555?135?mm?m2n?12n?12n?1?35?1llloln?2n?1?c1??1??n???3?2n?1??c2??0?n?5??0?2n?1?c??3???m??m??m?n?2n?1????2n?10???cn????ox2ox数值频散试验dxdz10mdt1ms10高阶差分为何会消除数值频散
常用滤波算法及C语言程序实现
A、方法：根据经验判断，确定两次采样允许的最大偏差值（设为A）每次检测到新值时判断：如果本次值与上次值之差<=A,则本次值有效如果本次值与上次值之差>A,则本次值无效,放弃本次值,用上次值代替本次值B、优点：能有效克服因偶然因素引起的脉冲干扰C、缺点无法抑制那种周期性的干扰平滑度差#define A 10 char value;char filter(){ char new_value; new_value = get_ad(); if ( ( new_value - value > A ) || ( value - new_value > A ) return value; return new_value; }2、中位值滤波法A、方法：连续采样N次（N取奇数）把N次采样值按大小排列取中间值为本次有效值B、优点：能有效克服因偶然因素引起的波动干扰对温度、液位的变化缓慢的被测参数有良好的滤波效果C、缺点：对流量、速度等快速变化的参数不宜#define N 11 char filter(){ char value_buf[N]; char count,i,j,temp; for ( count=0;count<N;count++) { value_buf[count] = get_ad(); delay(); } for (j=0;j<N-1;j++) { for (i=0;i<N-j;i++) { if ( value_buf>value_buf[i+1] ) { temp = value_buf; value_buf = value_buf[i+1]; value_buf[i+1] = temp; } } } return value_buf[(N-1)/2];}3、算术平均滤波法A、方法：连续取N个采样值进行算术平均运算N值较大时：信号平滑度较高，但灵敏度较低N值较小时：信号平滑度较低，但灵敏度较高N值的选取：一般流量，N=12；压力：N=4B、优点：适用于对一般具有随机干扰的信号进行滤波这样信号的特点是有一个平均值，信号在某一数值范围附近上下波动C、缺点：对于测量速度较慢或要求数据计算速度较快的实时控制不适用比较浪费RAM#define N 12 char filter(){ int sum = 0; for ( count=0;count<N;count++) { sum + = get_ad(); delay(); } return (char)(sum/N);}A、方法：把连续取N个采样值看成一个队列队列的长度固定为N每次采样到一个新数据放入队尾,并扔掉原来队首的一次数据.(先进先出原则)把队列中的N个数据进行算术平均运算,就可获得新的滤波结果N值的选取：流量，N=12；压力：N=4；液面，N=4~12；温度，N=1~4B、优点：对周期性干扰有良好的抑制作用，平滑度高适用于高频振荡的系统C、缺点：灵敏度低对偶然出现的脉冲性干扰的抑制作用较差不易消除由于脉冲干扰所引起的采样值偏差不适用于脉冲干扰比较严重的场合比较浪费RAM#define N 12 char value_buf[N];char i=0; char filter(){ char count; int sum=0; value_buf[i++] = get_ad(); if ( i == N ) i = 0; for ( count=0;count<N,count++) sum = value_buf[count]; return (char)(sum/N);}5、中位值平均滤波法（又称防脉冲干扰平均滤波法）A、方法：相当于“中位值滤波法”+“算术平均滤波法”连续采样N个数据，去掉一个最大值和一个最小值然后计算N-2个数据的算术平均值N值的选取：3~14B、优点：融合了两种滤波法的优点对于偶然出现的脉冲性干扰，可消除由于脉冲干扰所引起的采样值偏差C、缺点：测量速度较慢，和算术平均滤波法一样比较浪费RAM#define N 12 char filter(){ char count,i,j; char value_buf[N]; int sum=0; for (count=0;count<N;count++) { value_buf[count] = get_ad(); delay(); } for (j=0;j<N-1;j++) { for (i=0;i<N-j;i++) { if ( value_buf>value_buf[i+1] ) { temp = value_buf; value_buf = value_buf[i+1]; value_buf[i+1] = temp; } } } for(count=1;count<N-1;count++) sum += value[count]; return (char)(sum/(N-2));}6、限幅平均滤波法A、方法：相当于“限幅滤波法”+“递推平均滤波法”每次采样到的新数据先进行限幅处理，再送入队列进行递推平均滤波处理B、优点：融合了两种滤波法的优点对于偶然出现的脉冲性干扰，可消除由于脉冲干扰所引起的采样值偏差C、缺点：比较浪费RAM略参考子程序1、37、一阶滞后滤波法(低通滤波)A、方法：取a=0~1本次滤波结果=（1-a）*本次采样值+a*上次滤波结果B、优点：对周期性干扰具有良好的抑制作用适用于波动频率较高的场合C、缺点：相位滞后，灵敏度低滞后程度取决于a值大小不能消除滤波频率高于采样频率的1/2的干扰信号#define a 50 char value; char filter(){ char new_value; new_value = get_ad(); return (100-a)*value + a*new_value; }8、加权递推平均滤波法A、方法：是对递推平均滤波法的改进，即不同时刻的数据加以不同的权通常是，越接近现时刻的数据，权取得越大。
有限差分法地震正演模拟程序
有限差分法地震正演模拟程序一．二阶公式推导1．二维的弹性波动方程22222222x x x ZU U U U t x z x z λμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 22222222xz z z U U U U t z x x zλμμλμρρρ∂∂∂∂++=++∂∂∂∂∂ 2．对方程进行中间差分（1）首先对时间进行二阶差分()()()()()112,,,22222n n n xi k xi k xi kx x x x U U U U t t U t U t t U t t t +--++-+-∂==∂ ()()()()()112,,,22222n n n zi k zi k zi kz z z z U U U U t t U t U t t U t t t +--++-+-∂==∂ （2）对x 方向进行二阶差分()()()()()21,,1,22222n n n xi k xi k xi kx x x x U U U U x x U x U x x U x x x +--++-+-∂==∂ ()()()()()21,,1,22222n n n zi k zi k zi kz z z z U U U U x x U x U x x U x x x +--++-+-∂==∂ （3）对z 方向进行二阶差分()()()()()2,1,,122222n n nxi k xi k xi k x x x x U U U U z z U z U z z U z z z +--++-+-∂==∂ ()()()()()2,1,,122222n n nzi k zi k zi k z z z z U U U U z z U z U z z U z z z +--++-+-∂==∂ （4）对xz 进行差分2,1,11,11,11,11,111,11,11,11,122224n nzi k zi k z zn n n nzi k zi k zi k zi k n n n n zi k zi k zi k zi k U U U U x z x zx z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==2,1,11,11,11,11,111,11,11,11,122224n nxi k xi k x x n n n n xi k xi k xi k xi k n n n n xi k xi k xi k xi k U U U U x z x z x z U U U U U U U U x x z z x+-++-++--++-++---⎛⎫-∂∂∂∂⎛⎫== ⎪ ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭-----+==（5）把（1）（2）（3）（4）得到的结果带入波动方程3．写出差分方程()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n n n nxi k xi k xi kxi k xi k xi k nnnnnnnxi k xi k xi k zi k zi k zi k zi k U U U U U U t x U U UU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++()()()11,,,1,,1,22,1,,11,11,11,11,12222=2++4n n n nnnzi k zi k zi kzi k zi k zi knn n nnnnzi k zi k zi k xi k xi k xi k xi k U U U U U U t x U UUU UUUz xz λμρμλμρρ+-+-+-++-++----+-++-+--++即得到()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nxi k xi k xi kxi k xi k xi k n n n xi k xi k xi k n n n nzi k zi k zi k zi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭()()1,,1,,1,,122112,,,1,11,11,11,1222+=2-++4n n n n n nzi k zi k zi k zi k zi k zi k n n n zi k zi k zi k n n n nxi k xi k xi k xi k U U U U U U x z U U U t U U U U z x λμμρρλμρ+-+-+-++-++---⎛⎫-+-++ ⎪ ⎪ ⎪--++ ⎪ ⎪⎝⎭二．MATLAB 程序 clear; clc;Nx=200; Nz=300;fi=30;%%%主频t_step=300;%%%%时间采样点dx=10.0;%空间采样间隔——x 方向 dz=10.0;%空间采样间隔——z 方向 dt=0.001;%时间采样间隔——1mslambda=66*1e9;%砂岩拉梅常数lamdamu=44*1e9;%砂岩拉梅常数murho=2650;%砂岩密度%%%%%%模型扩展%%%%%vp=zeros(Nz+2,Nx+2);%纵波速度vs=zeros(Nz+2,Nx+2);%横波速度c=zeros(Nz+2,Nx+2);%交叉项系数for i=1:Nz+2for j=1:Nx+2vp(i,j)=sqrt((lambda+2*mu)/rho);vs(i,j)=sqrt((mu)/rho);c(i,j)=sqrt((lambda+mu)/rho);endend%%%%%% 震源%%%%%%%t_wavelet=[1:t_step]*dt-1.0/fi;source=((1-2*(pi*fi*t_wavelet).^2).*exp(-(pi*fi*t_wavelet).^2));% 雷克子波amp=sqrt(2.0);% 振幅source_x=floor(Nx/2)+1;% 震源位置——x坐标source_z=2;% 震源位置——z坐标source_amp=zeros(Nz+2,Nx+2);% 震源振幅初始化（所有点处均为0）source_amp(source_z,source_x)=amp;% 震源振幅，在位置source_z,source_x处振幅为amp，其它位置为0 %%%%%%%%%%%%%%%%%%%%%%%%%%%U=zeros(Nz,Nx);% 弹性波x分量V=zeros(Nz,Nx);% 弹性波z分量U0=zeros(Nz+2,Nx+2);% 前一时刻的UU1=zeros(Nz+2,Nx+2);% 当前时刻的UU2=zeros(Nz+2,Nx+2);% 下一时刻的UV0=zeros(Nz+2,Nx+2);% 前一时刻的VV1=zeros(Nz+2,Nx+2);% 当前时刻的VV2=zeros(Nz+2,Nx+2);% 下一时刻的Vrecord_u=zeros(t_step,Nx);% x方向接收到的地震记录——Urecord_v=zeros(t_step,Nx);% x方向接收到的地震记录——V%%%%%% 有限差分计算U V %%%%%%for it=1:t_stepfor x=2:Nx+1for z=2:Nz+1U2(z,x)=2*U1(z,x)-U0(z,x)+(dt*dt)*(vp(z,x)^2*(1.0/(dx*dx))*(U1(z,x+1)-2*U1(z,x) +U1(z,x-1))+vs(z,x)^2*(1.0/(dz*dz))*(U1(z+1,x)-2*U1(z,x)+U1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(V1(z+1,x+1)-V1(z+1,x-1)-V1(z-1,x+1)+V1(z-1,x-1)))+source(it)*sour ce_amp(z,x)*dt*dt;V2(z,x)=2*V1(z,x)-V0(z,x)+(dt*dt)*(vs(z,x)^2*(1.0/(dx*dx))*(V1(z,x+1)-2*V1(z,x) +V1(z,x-1))+vp(z,x)^2*(1.0/(dz*dz))*(V1(z+1,x)-2*V1(z,x)+V1(z-1,x))+c(z,x)^2*(1. 0/(4*dx*dz))*(U1(z+1,x+1)-U1(z+1,x-1)-U1(z-1,x+1)+U1(z-1,x-1)));endendU0=U1;U1=U2;V0=V1;V1=V2;record_u(it,:)=U2(2,[2:201]);record_v(it,:)=V2(2,[2:201]);U=U2(2:301,2:201);V=V2(2:301,2:201);end%%%%% 绘制U 和V 的波场图%%%%%%figureimagesc(U(1:300,1:200));title('U');xlabel('x');ylabel('z');figureimagesc(V(1:300,1:200));title('V');xlabel('x');ylabel('z');%%%%% 绘制U 和V 的地震记录%%%%%%figureimagesc(record_u(1:300,1:200));title('U');xlabel('x');ylabel('t_step');figureimagesc(record_v(1:300,1:200));title('V');xlabel('x');ylabel('t_step');U分量波场图V分量波场图U分量地震记录V分量地震记录。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
#include<stdio.h>
#include<math.h>
#define fm 30
#define dt 0.001
#define PI 3.1415926
#define Nt 401
#define Nx 200
#define Nz 200
//---------------加载震源，雷克子波-----------------------------
void fun(float source[])
{
FILE *fp;
int it,i;
float t1,t2,t0;
for(i=0;i<Nt;i++)
source[i]=0.0;
t0 = 1.0/fm;
printf("t0 = %f\n", t0);
fp=fopen("ricker.txt","w");
for(it=0;it<Nt;it++)
{
t1=it*dt;
t2=t1-t0;
source[it]=(1.0-2.0*pow(PI*fm*t1,2.0))*exp(-pow(PI*fm*t1,2.0));
fprintf(fp,"%.8f %.8f\n",t1,source[it]);
}
fclose(fp);
}
//--------------------主程序，差分算子----------------------------------------
void main()
{ FILE *fp, *fpwf;
int it,ix,iz, order,N;
int nsx,nsz;
float a0,a1,a2,a3,a4,a5,a6;
int flag;
float dx,dz;
float Un[Nx][Nz], Unm[Nx][Nz], Unp[Nx][Nz], v[Nx][Nz],record[Nx][Nt];//Un表示n时刻，Unm表示n-1时刻，Unp表示n+1时刻波场
float source[Nt];
nsx=Nx/2;
nsz=0.0;//震源位置
N=6;
order=2*N;//精度
if(order == 2)
{
a0=-2.00000;a1=1.00000;a2=0;a3=0;a4=0;a5=0;a6=0;
}
if(order == 4)
{
a0=-2.50000;a1=1.33333;a2=-8.33333e-2;a3=0;a4=0;a5=0;a6=0;
}
if(order == 6)
{
a0=-2.72222;a1=1.50000;a2=-1.50000e-1;a3=1.11111e-2;a4=0;a5=0;a6=0;
}
if(order == 8)
{
a0=-2.84722;a1=1.60000;a2=-2.00000e-1;a3=2.53968e-2;a4=-1.78571e-3;a5=0;a6=0;
}
if(order == 10)
{
a0=-2.92722;a1=1.66667;a2=-2.38095e-1;a3=3.96825e-2;a4=-4.96032e-3;a5=3.17460e-4;a6=0;
}
if(order == 12)
{
a0=-2.98278;a1=1.71429;a2=-2.67857e-1;a3=5.29101e-2;a4=-8.92857e-3;a5=1.03896e-3;a6=-6.01251e-5;
}
/*for (ix=0;ix<Nx;ix++)
{ for (iz=0;iz<(Nz*2/5);iz++)
v[ix][iz]=2000;
for (iz=(Nz*2/5);iz<(Nz*3/5);iz++)
v[ix][iz]=3000;
for (iz=(Nz*3/5);iz<Nz;iz++)
v[ix][iz]=4000;//给速度赋值
}
*/
for (ix=0;ix<Nx;ix++)
{ for (iz=0;iz<Nz/2;iz++)
v[ix][iz]=4000;
for (iz=Nz/2;iz<Nz;iz++)
v[ix][iz]=4000;
}
dx=10.0;
dz=10.0;//赋值
for(ix=0;ix<Nx;ix++)
{
for(iz=0;iz<Nz;iz++)
{
Un[ix][iz]=0.0;
Unm[ix][iz]=0.0;
Unp[ix][iz]=0.0;
}
}
fun(source);
fpwf = fopen("Wavefieldall.dat", "wb");
for(it=0;it<Nt;it++)
{
for(ix=6;ix<(Nx-6);ix++)
for(iz=6;iz<(Nz-6);iz++)
{
Unp[ix][iz]=2*Un[ix][iz]-
Unm[ix][iz]+0.5*v[ix][iz]*v[ix][iz]*dt*dt/(dx*dx)*(a0*Un[ix][iz]+a1*(Un[ix+1][iz]+Un[ix-1][iz])+a2*(Un[ix+2][iz]+Un[ix-2][iz])+a3*(Un[ix+3][iz]+Un[ix-3][iz])+a4*(Un[ix+4][iz]+Un[ix-4][iz])+a5*(Un[ix+5][iz]+Un[ix-5][iz])+a6*(Un[ix+6][iz]+Un[ix-6][iz]))
+0.5*v[ix][iz]*v[ix][iz]*dt*dt/(dz*dz)*(a0*Un[ix][iz]+a1*(Un[ix][iz+1]+Un[ix][iz-
1])+a2*(Un[ix][iz+2]+Un[ix][iz-2])+a3*(Un[ix][iz+3]+Un[ix][iz-3])+a4*(Un[ix][iz+4]+Un[ix][iz-4])+a5*(Un[ix][iz+5]+Un[ix][iz-5])+a6*(Un[ix][iz+6]+Un[ix][iz-6]));
}
Unp[nsx][nsz]+=source[it];
for(ix=0;ix<Nx;ix++)
{
for(iz=0;iz<Nz;iz++)
{Unm[ix][iz]=Un[ix][iz];
Un[ix][iz]=Unp[ix][iz];}
}
for(ix=0;ix<Nx;ix++)
{
record[ix][it]=Un[ix][nsz];
}
//------------------------记录文件-------------------------
//-------记录Un文件-----------
for(ix=0; ix<Nx; ix++)
{
fwrite(Un[ix], sizeof(float), Nz, fpwf);
}
if(it==400)
{
fp=fopen("wavefield.dat","wb");
for(ix=0; ix<Nx; ix++)
{
fwrite(Un[ix], sizeof(float), Nz, fp);
}
fclose(fp);
}
if(it%20==0) printf("it=%d\n", it);
}
fclose(fpwf);
//--------记录record文件-------
fp=fopen("record.dat","wb");
for(ix=0;ix<Nx;ix++)
{
fwrite(record[ix],sizeof(float),Nt,fp);
}
fclose(fp);
fp=fopen("wave1.txt","w");
for(it=0;it<Nt;it++)
{
fprintf(fp,"%.8f\n",Unp[nsx][it]);
}
fclose(fp);
}。