基于MATLAB的地震数据的分析

基于MATLAB的地震数据的分析

孙玉柱冯光房桂梅

摘要：地震波原始数据中存在的干扰信号，会影响震相分析的准确性。为了滤除干扰信号，对地震波原始信号进行了频谱分析，给出了一种基于MATLAB的FIR数字滤波器的优化设计方案，将其用于地震波数据的分析中，并进行了仿真分析。仿真结果表明，FIR数字滤波器对地震波原始信号进行滤波处理后，提高了震相分析的准确性，得到了理想的效果，达到了预期的目的。

关键词：MATLAB；FIR数字滤波器；优化；滤波

the Analysis of Earthquake Data Based on MATLAB

SUN Yuzhu,FENG Guang,FANG Guimei

Abstract: The interference that existed in the earthquake data will affect the accuracy of the seismic phase analysis. In order to filter the disturbance signal, this paper carries out spectrum analysis of the earthquake data, proposes an optimum design method for FIR digital filter based on MATLAB and applies it to the analysis of earthquake data. After the filter of the noise jamming, the true information of the earthquake wave is clearly reflected. The simulation results manifest that it can

improve the accuracy of seismic phase analysis and arrive at the purpose desired.

Key words: MATLAB；FIR digital filter；optimization；filter

1 引言

地震带给人类的损失是巨大的，汶川大地震依旧在我们的记忆深处清晰存在。大地震的每次不约而至，都对国家和人民造成了巨大的损失。地震预测是世界性难题，全世界的地震科学家不断在探索，尽最大努力减少其破坏性。地震台站提供的地震观测资料的可靠性和准确性，是地震学家进行地震预测的基础[1]。但是地震波信号变化的不平稳性和复杂性给地震的分析和预测带来了很大困难，并且在地震波的原始记录中往往还掺杂着来自外界的各种干扰，如仪器、环境噪声、爆破、采矿、火车的震动等，这些都给地震波的分析带来了严重影响，甚至导致分析结果的错误。为保证地震分析的准确性，可先对原始记录进行频谱分析，选择性能优良的滤波器对其进行优化处理，把干扰信号尽量滤除，然后再对处理后的地震波数据进行分析处理，则会得到良好的效果。

地震记录的数字化使得利用计算机对地震信号进行分析处理得以实现。在数字信号的分析处理中，Fourier变换和数字滤波器的应用极为普遍，语音、雷达、地震、图像、机械振动、地质勘探等众多领域都广泛采用数字滤波器。MATLAB是一种集数值分析、矩阵运算、信号处理和图像显示于一体、功能极其强大的高性能软件，其工具箱

中包含了各种经典和现代数字信号处理技术，很容易实现Fourier 变换和各种数字滤波器的设计，在地震数据的分析处理中起着重要作用。

本文首先介绍了数字滤波器，给出了基于MATLAB 的FIR 数字滤波器的一种优化设计方案，最后将其用于地震波数据的处理中，并进行了仿真分析。

2快速Fourier 变换（FFT ）及频谱分析

在地震波的原始记录数据中往往夹杂不同频率范围的噪声干扰信号，为了显示出地震波数据中的优势频率和干扰频率，保证地震分析的准确性，应首先采用频谱分析，再针对干扰波的频率范围，设计合适的滤波器参数。

在对有限长信号序列进行频谱分析时，离散Fourier 变换（DFT ）应用非常广泛，它可以很好地反映序列的频谱特性[2]。

设()x n 是一个长度为N 的有限长序列，则()x n 的N 点离散傅里叶变换为[3] ()1

20()(), 0,1,...,1N j N kn n X k x n e k N π--===-∑ （1）

DFT 是信号分析与处理中的一种重要变换，但是当N 较大时，其计算量太大。快速Fourier 变换（FFT ）是减少 DFT 运算次数的一种快速算法，通过在时域将序列逐次分解为一组子序列，然后利用子序列的DFT 来实现整个序列的DFT ，从而减少离散Fourier 变换的运算量，提高了计算效率。

在MATLAB 中可调用函数()f t X =fft x ,N 来进行快速傅里叶变换[4]。

其中，

x为时域内的输入信号序列，N为序列长度，f X为频率域的输

t

出信号，即

x的频谱特征。

t

3 FIR数字滤波器的优化设计

3.1 数字滤波器的选择

在地震分析中必须要先对原始信号做滤波处理，滤波的目的是为了去除噪声，使原始信号通过滤波器后能够清晰地显示出优势频率，为更好的分析地震信号（比如震相等）做准备。滤波器包括模拟滤波器和数字滤波器，模拟滤波器又可分为无源滤波器（主要由R、C和L构成）和有源滤波器（主要由集成运放和R、C元件构成），数字滤波器可用计算机软件或大规模集成数字硬件实现。模拟滤波器存在电压漂移、温度漂移和噪声等问题，而数字滤波器不存在这些问题，可以达到很高的稳定度和精度。

根据实现的网络结构不同，数字滤波器可分为无限脉冲响应（IIR）滤波器和有限脉冲响应（FIR）滤波器。

1. IIR数字滤波器

IIR数字滤波器最大的优点是可取得非常好的通带和阻带衰减，还可得到准确的通带与阻带的边缘频率，而且滤波时需要的计算量小；缺点是滤波器的响应灵活性较差，不具有线性相位（若得到线性相位则需要用全通系统进行相位补偿）且存在稳定性问题，而这些在实际应用（如图像信号处理、地震信号处理、数据通信等领域）中又是非常重要的。

2. FIR数字滤波器

FIR 数字滤波器又称为卷积滤波器，通常采用迭代算法来满足设定的技术指标，由于FIR 系统是全零点系统，其单位抽样响应h(n)为有限长，因此容易实现某种对称性，从而获得在设计任意幅频特性的同时保证严格的线性相位特性；由于采用非递归结构，则不存在稳定性问题，运算误差也较少，并且能容易实现多通带、多阻带的滤波器的设计；又因为FIR 数字滤波器在通带内具有恒定的幅频特性和线性相位特性，从而使信号通过时不失真，有时还可实现零相位滤波。

考虑到地震波数据的特点，本文选用FIR 数字滤波器。

3.2 FIR 数字滤波器的MATLAB 优化设计

FIR 数字滤波器的传递函数为[5] 10

()()()()N n n Y z H z h n z X z --===∑ （2） 其中，()h n 是FIR 滤波器的单位脉冲响应，其长度为N ，非零区间为

[0,1]N -。（1）式表明()H z 是1z -的(1)N -次多项式，它在z 平面上有(1)N -个零点。原点0z =是(1)N -阶重极点。因此，()H z 肯定是稳定的。 由FIR 滤波器的传递函数，可得到其频率响应表达式

1

0()()N j j n n H e h n e ω

ω--==∑ （3） 如果()h n 是实序列，并且满足下列中心对称条件[6]

()(1)h n h N n =--或()(1)h n h N n =--- （4） 则FIR 滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。

FIR 滤波器虽然具有相位滞后的缺点，但是其相位滞后和群延迟

在整个频带上是相等且不变的。一个N 阶的线性相位FIR 滤波器群延迟为2N ，即滤波后的信号简单地延长常数个时间步长，这一特性使通带频率内信号通过滤波器后仍然保持原有波形形状而无相位失真。

FIR 数字滤波器的设计方法很多，本文采用目前最优秀的设计方法——切比雪夫逼近法。切比雪夫逼近法是一种等波纹逼近法，它使误差在整个频带均匀分布，对同样的技术指标，这种逼近法需要的滤波器阶数低，而对同样的滤波器阶数，这种逼近法的最大误差最小。设希望设计的滤波器幅度特性为()d H ω，实际设计的滤波器幅度特性为()g H ω，误差加权系数为()W ω，则切比雪夫最佳一致逼近准则是使加权函数[7]

()()()()d g E W H H ωωωω??=-?? （5） 的最大值达到最小，即满足 min max ()A E ωω∈????

（6） 切比雪夫逼近理论可以解决()E ω的存在性、唯一性及如何构造等一系列问题，克服了通带和阻带的边缘不易精确确定和Gibbs 现象等多种缺点，使设计出的滤波器具有明显优点。因为在一定意义上对()d H ω做最佳逼近，可获取较好的通带和阻带性能，并能准确地指定通带和阻带的边缘频率，此时又具有了IIR 滤波器的优点，所以这是一种很好的设计方法。

MATLAB 提供了大量设计FIR 数字滤波器的函数，下面针对切比雪夫逼近法，介绍其设计步骤和方法。

（1）根据得出的地震波原始记录的频谱分析波形，确定FIR 数

字滤波器的技术指标；

（2）用函数[]()

M,F,A,W=remezord f,a,dev,Fs[8]估算设计的等波纹逼

00

近法的参数：最低滤波器阶数M、频率向量

F、幅度向量0A和加权

向量W。其中，Fs为采样频率；f可以是模拟频率或归一化频率，但必须以0开始，以Fs2（用归一化频率为1）结束，而且省略了0和Fs2两个频点；dev为各逼近段允许的幅度响应偏差（波纹振幅）；a为滤波器在各个频段上的幅度值，一般对通带取值为1，对阻带取值为0。

（3）用函数()

h=remez M,F,A,W[8]完成等纹波FIR滤波器的设计，

00

并用函数()

y=filtfilt h,1,x完成对输入信号x的滤波，得到y=filter h,1,x或()

输出信号y。

4 基于MATLAB的地震数据的分析

下面以汕头台受到距台站300m处的汽车干扰的波形记录图的数字地震记录资料为例，进行基于MATLAB的地震数据的分析。该地震台采用的是频带范围为0.05~20Hz 的FBS-3A型宽频带数字地震记录仪和EDAS-C24型数据采集器，系统的采样频率为50Hz。

将该地震波原始数据保存在MATLAB的work文件夹中，并对其进行地频谱分析，如图1所示。从图1中可以很清楚的看出地震波原始记录中地震信号的优势频率为0.25Hz，频段范围在0~1Hz之间；干扰的优势频率为12.5Hz，频段范围在10~15Hz之间。

50100150

-10000

1000时间/s 振幅/c o u n t s 时域波形图

0510

152025

012x 10

6频率/Hz

频谱密度频域波形图0510

152025100

10510

10频率/Hz 频谱密度取对数频域波形图

图1 地震波原始数据的时域波形图和频谱分析图

Figure 1 the original time-domain and spectrum analysis of the

earthquake data

为了滤除干扰信号，最大限度的保留其中的有用信号，选取的FIR 滤波器为带阻滤波器，其参数选定为：通带上截止频率Fp1=7，阻带下截止频率Fs1=7.1，阻带上截止频率Fs1=18.9，通带下截止频率Fp2=19，通带波纹峰值dp=0.01，阻带波纹峰值ds=0.01。则设计的带阻FIR 滤波器的频率特性曲线如图2所示。

051015

2025-6-4

-20

4Frequency (Hz)P h a s e (d e g r e e s )0

510152025

-150-100

-50

50

Frequency (Hz)M a g n i t u d e (d B )

图2 带阻FIR 滤波器的频率特性曲线

Figure 2 the frequency characteristic of band-stop FIR filter

由图2可以看出设计的FIR 滤波器的特性满足设计要求。利用该滤波器对地震波原始记录信号进行滤波处理，滤波前后的时域波形图如图3所示，图4为滤波前后频域图的对比。

050

100150

时间/s

振幅/c o u n t s 滤波前时域波形图

050

100150

时间/s 振幅/c o u n t s 滤波后时域波形图

图3 滤波前后地震波时域波形图的比较

Figure 3 the comparison of time-domain between before and after filter

processing

51015202505

1015x 10

5频率/Hz 频谱密度滤波前频域波形图

0510

1520250510

15

5频率/Hz 频谱密度

滤波后频域波形图

图4 滤波前后频域波形图的比较

Figure 4 the comparison of frequency-domain between before and after

filter processing

由图3和图4可以看出，地震波原始记录信号经FIR 滤波器处理后，干扰波信号被滤除了，而地震波信号很好的显示了出来。

由于FIR 滤波器具有相位滞后的特点，所以滤波后的信号和原信号相比有了一定程度的相位延迟。为了不产生相位延迟，使滤波后的波形图能更好的与原始信号进行比较，可以对程序2（见附件）进行修改：

第38行后加 m=(M-1)/2/Fs; %计算相位延迟 第50行后加 t=t-m; %计算相位延迟 第53行后加 xlim([0,max(t)]);

然后再对地震波原始记录信号进行处理，则滤波前后地震波时域波形图如图5所示，该图形非常明显的反映出了滤波前后地震波的特

征，得到了理想的效果，达到了预期设计的目的。

050

100150

-1000-500

500

时间/s

振幅/c o u n t s 滤波前时域波形图

0204060

80100120140-600-400

-2000

200

时间/s 振幅/c o u n t s 滤波后时域波形图

图5滤波前后时域波形图的比较

Figure 5 the comparison of time-domain between before and after filter

processing

5 结束语

地震数据中包含了很多干扰信号，直接影响分析的准确性。干扰信号的处理虽然已经有很多方法，但本文采用目前非常流行的MATLAB 软件，使用数字信号处理中的快速Fourier 变换和最优滤波器的设计方法，对采集到的地震数据进行频谱分析和滤波处理，去除干扰并最大限度地保留有用的频率成份，并使信号无失真，提高了震相分析的准确度。该方法可用于结构地震动力分析、地震台等领域，对地震的观测、分析、预报和研究有着极其重要的作用。

参考文献

[1] 宋建锁. 滤波在地震分析中的应用[J]. 防灾技术高等专科学校学报, 2006,8(1):75-79

[2]李敬,甘延锋,黄友明. 数字地震记录中干扰波的排除[J]. 防灾技术高等专科学校学报, 2004,6(3):20-25

[3] 丁玉美, 高西金. 数字信号处理[M]. 西安:西安电子科技大学出版社,2004

[4] 万永革. 数字信号处理的MATLAB实现[M]. 北京:科学出版社, 2007

[5]胡广书. 数字信号处理[M]. 北京:清华大学出版社, 2003

[6]陈鲁刚,王锐锋. 用MATLAB实现滤波器对数字地震资料处理初探[J]. 防灾技术高等专科学校学报,2006,8(1):62-65

[7] 陈后金. 数字信号处理[M]. 北京:高等教育出版社, 2004

[8] 唐向红, 岳恒立, 郑雪峰. MATLAB及在电子信息类课程中的应用[M]. 北京:电子工业出版社, 2006

附件

程序1 对地震记录数据进行频谱分析的实现程序

load i.txt %调入原始数据文件（加载地震波数据）Xt=i; %得到原始信号序列

Fs=50; %采样频率为50Hz

dt=1/Fs; %采样间隔（单位为s）

N=length(Xt); %原始信号序列长度

Xf=fft(Xt); %对原始波形数据进行快速Fourier变换

n=0:N-1;

t=n*dt; %得到时间序列

f=n/(N*dt); %得到频率序列

subplot(3,1,1); %时域坐标方框图

plot(t, Xt); %画出时域中的原始波形图xlabel('时间/s'); %X轴标示

ylabel('振幅/counts'); %Y轴标示

title('时域波形图'); %加注标题

grid on

subplot(3,1,2); %频域坐标方框图

plot(f,abs(Xf)); %画出频域中的FFT波形图

xlabel('频率/Hz'); %X轴标示

ylabel('频谱密度'); % Y轴标示

title('频域波形图'); %加注标题

xlim([0 Fs/2]); %频域之画出采样频率的一半grid on

subplot(3,1,3); %频域坐标方框图semilogy(f,abs(Xf)); %画出频域中的FFT波形图，Y轴为对数xlabel('频率/Hz'); %X轴标示

ylabel('频谱密度取对数'); % Y轴标示

title('频域波形图'); %加注标题

xlim([0 Fs/2]); %频域之画出采样频率的一半grid on

程序2 用FIR数字滤波器实现的滤波程序

load i.txt %调入原始数据文件（加载地震波数据）Xt=i; %得到原始信号序列

Fs=50; %采样频率为50Hz

dt=1/Fs; %采样间隔（单位为s）

N=length(Xt); %原始信号序列长度

Xf=fft(Xt); %对原始波形数据进行快速Fourier变换n=0:N-1;

t=n*dt; %得到时间序列

f=n/(N*dt); %得到频率序列

figure(1)

subplot(2,1,1); %时域坐标方框图

plot(t, Xt); %画出时域中的原始波形图

xlabel('时间/s'); %X轴标示

ylabel('振幅/counts'); %Y轴标示

title('时域波形图'); %加注标题

grid on

subplot(2,1,2); %频域坐标方框图

plot(f,abs(Xf)); %画出频域中的FFT波形图xlabel('频率/Hz'); %X轴标示

ylabel('频谱密度'); % Y轴标示

title('频域波形图'); %加注标题

xlim([0 Fs/2]); %频域之画出采样频率的一半grid on

Fp1=7; %带阻滤波器通带上截止频率Fs1=7.1; %带阻滤波器阻带下截止频率Fs2=18.9; %带阻滤波器阻带上截止频率Fp2=19; %带阻滤波器通带下截止频率dp=0.01; %带阻滤波器通带波纹峰值

ds=0.01; %带阻滤波器阻带波纹峰值W1=2*Fp1/Fs;

W2=2*Fs1/Fs;

W3=2*Fs2/Fs;

W4=2*Fp2/Fs; %通带和阻带归一化频率f1=[W1 W2 W3 W4]; %频率向量

a=[1 0 1]; %振幅向量

dev=[dp ds dp]; %各逼近段允许的幅度响应偏差[M,Fo,Ao,W]=remezord(f1,a,dev)；%返回设计FIR带阻滤波器的参数

h=remez(M,Fo,Ao,W); %设计FIR带阻滤波器figure(2);

freqz(h,1,N,Fs); %求出滤波器幅频相频特性grid on;

figure(3);

subplot(2,1,1);

plot(t,Xt); %画出时域中的原始波形图xlabel('时间/s'); %X轴标示

ylabel('振幅/counts'); %Y轴标示

title('滤波前时域波形图'); %加注标题

grid on

z=filter(h,1,Xt); %对原始信号进行滤波subplot(2,1,2);

plot(t,z); %画出滤波后的时域波形图xlabel('时间/s'); %X轴标示

ylabel('振幅/counts'); %Y轴标示

title('滤波后时域波形图'); %加注标题

grid on

Q=length(z); %求滤波后信号序列长度figure(4);

subplot(2,1,1);

plot(f,abs(Xf)); %画出频域中的原始FFT波形图xlabel('频率/Hz'); %X轴标示

ylabel('频谱密度'); % Y轴标示

title('滤波前频域波形图'); %加注标题

xlim([0 Fs/2]); %频域直画出采样频率的一半grid on

p=fft(z); %对滤波后的波形数据进行快速Fourier变换subplot(2,1,2);

plot([0:Q-1]/(Q*dt),abs(p)); %画出频域中的滤波后信号的FFT 波形图

xlabel('频率/Hz'); %X轴标示

ylabel('频谱密度'); % Y轴标示

title('滤波后频域波形图'); %加注标题

xlim([0 Fs/2]); %频域直画出采样频率的一半

grid on

matlab在统计数据的描述性分析的应用

统计数据的描述性分析 一、实验目的 熟悉在matlab中实现数据的统计描述方法，掌握基本统计命令：样本均值、样本中位数、样本标准差、样本方差、概率密度函数pdf、概率分布函数df、随机数生成rnd。 二、实验内容 1 、频数表和直方图 数据输入,将你班的任意科目考试成绩输入 >> data=[91 78 90 88 76 81 77 74]; >> [N,X]=hist(data,5) N = 3 1 1 0 3 X = 75.7000 79.1000 82.5000 85.9000 89.3000 >> hist(data,5)

2、基本统计量 1) 样本均值 语法: m=mean(x) 若x 为向量，返回结果m是x 中元素的均值； 若x 为矩阵，返回结果m是行向量，它包含x 每列数据的均值。 2) 样本中位数 语法: m=median(x) 若x 为向量，返回结果m是x 中元素的中位数； 若x 为矩阵，返回结果m是行向量，它包含x 每列数据的中位数3) 样本标准差 语法:y=std(x) 若x 为向量，返回结果y 是x 中元素的标准差； 若x 为矩阵，返回结果y 是行向量，它包含x 每列数据的标准差

std(x)运用n-1 进行标准化处理，n是样本的个数。 4) 样本方差 语法：y=var(x); y=var(x,1) 若x 为向量，返回结果y 是x 中元素的方差； 若x 为矩阵，返回结果y 是行向量，它包含x 每列数据的方差 var(x)运用n-1 进行标准化处理（满足无偏估计的要求），n 是样本的个数。var(x,1)运用n 进行标准化处理，生成关于样本均值的二阶矩。 5) 样本的极差（最大之和最小值之差） 语法：z= range(x) 返回结果z是数组x 的极差。 6) 样本的偏度 语法：s=skewness(x) 说明：偏度反映分布的对称性，s>0 称为右偏态，此时数据位于均值右边的比左边的多；s<0,情况相反；s 接近0 则可认为分布是对称的。 7) 样本的峰度 语法：k= kurtosis(x) 说明：正态分布峰度是3，若k 比3 大得多，表示分布有沉重的尾巴，即样本中含有较多远离均值的数据，峰度可以作衡量偏离正态分布的尺度之一。 >> mean(data) ,

MATLAB数据分析与多项式计算(M)

第7章 MATLAB数据分析与多项式计算 6.1 数据统计处理 6.2 数据插值 6.3 曲线拟合 6.4 离散傅立叶变换 6.5 多项式计算 6.1 数据统计处理 6.1.1 最大值和最小值 MATLAB提供的求数据序列的最大值和最小值的函数分别为max 和min，两个函数的调用格式和操作过程类似。 1．求向量的最大值和最小值 求一个向量X的最大值的函数有两种调用格式，分别是： (1) y=max(X)：返回向量X的最大值存入y，如果X中包含复数元素，则按模取最大值。 (2) [y,I]=max(X)：返回向量X的最大值存入y，最大值的序号存入I，如果X中包含复数元素，则按模取最大值。 求向量X的最小值的函数是min(X)，用法和max(X)完全相同。 例6-1 求向量x的最大值。 命令如下： x=[-43,72,9,16,23,47]; y=max(x) %求向量x中的最大值 [y,l]=max(x) %求向量x中的最大值及其该元素的位置 2．求矩阵的最大值和最小值 求矩阵A的最大值的函数有3种调用格式，分别是： (1) max(A)：返回一个行向量，向量的第i个元素是矩阵A的第i 列上的最大值。 (2) [Y,U]=max(A)：返回行向量Y和U，Y向量记录A的每列的最大值，U向量记录每列最大值的行号。 (3) max(A,[],dim)：dim取1或2。dim取1时，该函数和max(A)完全相同；dim取2时，该函数返回一个列向量，其第i个元素是A矩阵的第i行上的最大值。 求最小值的函数是min，其用法和max完全相同。

例6-2 分别求3×4矩阵x中各列和各行元素中的最大值，并求整个矩阵的最大值和最小值。 3．两个向量或矩阵对应元素的比较 函数max和min还能对两个同型的向量或矩阵进行比较，调用格式为： (1) U=max(A,B)：A,B是两个同型的向量或矩阵，结果U是与A,B 同型的向量或矩阵，U的每个元素等于A,B对应元素的较大者。 (2) U=max(A,n)：n是一个标量，结果U是与A同型的向量或矩阵，U的每个元素等于A对应元素和n中的较大者。 min函数的用法和max完全相同。 例6-3 求两个2×3矩阵x, y所有同一位置上的较大元素构成的新矩阵p。 6.1.2 求和与求积 数据序列求和与求积的函数是sum和prod，其使用方法类似。设X是一个向量，A是一个矩阵，函数的调用格式为： sum(X)：返回向量X各元素的和。 prod(X)：返回向量X各元素的乘积。 sum(A)：返回一个行向量，其第i个元素是A的第i列的元素和。 prod(A)：返回一个行向量，其第i个元素是A的第i列的元素乘积。 sum(A,dim)：当dim为1时，该函数等同于sum(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素之和。 prod(A,dim)：当dim为1时，该函数等同于prod(A)；当dim为2时，返回一个列向量，其第i个元素是A的第i行的各元素乘积。 例6-4 求矩阵A的每行元素的乘积和全部元素的乘积。 6.1.3 平均值和中值 求数据序列平均值的函数是mean，求数据序列中值的函数是median。两个函数的调用格式为： mean(X)：返回向量X的算术平均值。 median(X)：返回向量X的中值。

Matlab对采样数据进行频谱分析

使用Matlab对采样数据进行频谱分析 1、采样数据导入Matlab 采样数据的导入至少有三种方法。 第一就是手动将数据整理成Matlab支持的格式，这种方法仅适用于数据量比较小的采样。 第二种方法是使用Matlab的可视化交互操作，具体操作步骤为：File --> Import Data，然后在弹出的对话框中找到保存采样数据的文件，根据提示一步一步即可将数据导入。这种方法适合于数据量较大，但又不是太大的数据。据本人经验，当数据大于15万对之后，读入速度就会显著变慢，出现假死而失败。 第三种方法，使用文件读入命令。数据文件读入命令有textread、fscanf、load 等，如果采样数据保存在txt文件中，则推荐使用 textread命令。如 [a,b]=textread('data.txt','%f%*f%f'); 这条命令将data.txt中保存的数据三个三个分组，将每组的第一个数据送给列向量a，第三个数送给列向量b，第二个数据丢弃。命令类似于C语言，详细可查看其帮助文件。文件读入命令录入采样数据可以处理任意大小的数据量，且录入速度相当快，一百多万的数据不到20秒即可录入。强烈推荐！ 2、对采样数据进行频谱分析 频谱分析自然要使用快速傅里叶变换FFT了，对应的命令即 fft ，简单使用方法为：Y=fft(b,N)，其中b即是采样数据，N为fft数据采样个数。一般不指定N，即简化为Y=fft(b)。Y即为FFT变换后得到的结果，与b的元素数相等，为复数。以频率为横坐标，Y数组每个元素的幅值为纵坐标，画图即得数据b的幅频特性；以频率为横坐标，Y数组每个元素的角度为纵坐标，画图即得数据b的相频特性。典型频谱分析M程序举例如下： clc fs=100; t=[0:1/fs:100]; N=length(t)-1;%减1使N为偶数 %频率分辨率F=1/t=fs/N p=1.3*sin(0.48*2*pi*t)+2.1*sin(0.52*2*pi*t)+1.1*sin(0.53*2*pi*t)... +0.5*sin(1.8*2*pi*t)+0.9*sin(2.2*2*pi*t); %上面模拟对信号进行采样，得到采样数据p，下面对p进行频谱分析 figure(1) plot(t,p); grid on title('信号 p(t)'); xlabel('t') ylabel('p')

Matlab大数据处理

Matlab大数据处理2：硬盘访问.mat文件 分类：Matlab Hack2013-09-08 20:16 146人阅读评论(0) 收藏举报Matlab程序中经常要访问.mat文件，通常在作法是用load函数直接加载.mat文件。如果.mat文件非常大，超过了系统可用内存的时候该怎么办呢？Matlab2013b为提供了matfile函数，matfile函数可以通过索引直接访问.mat文件中的Matlab变量，而无需将.mat文件加载入内存。 matfile有两种用法： m = matfile(filename)，用文件名创建matfile对象，通过这个对象可以直接访问mat文件中的matlab变量。 m = matfile(filename,'Writable',isWritable)，isWritable开启或关闭文件写操作。 使用示例： 1. 向mat文件中写入变量 x = magic(20); m = matfile('myFile.mat'); % 创建一个指向myFile.mat的matfile对象 m.x = x; % 写入x m.y(81:100,81:100) = magic(20); % 使用坐标索引

2. 加载变量 filename = 'topography.mat'; m = matfile(filename); topo = m.topo; %读取变量topo [nrows,ncols] = size(m,'stocks'); %读取stocks变量的size avgs = zeros(1,ncols); for idx = 1:ncols avgs(idx) = mean(m.stocks(:,idx)); end 3. 开启写权限 filename = 'myFile.mat'; m = matfile(filename,'Writable',true); 或者 m.Properties.Writable = true;

第6章matlab数据分析与多项式计算_习题答案

第6章 MATLAB数据分析与多项式计算 习题6 一、选择题 1．设A=[1,2,3,4,5;3,4,5,6,7]，则min(max(A))的值是（）。B A．1 B．3 C．5 D．7 2．已知a为3×3矩阵，则运行mean(a)命令是（）。B A．计算a每行的平均值 B．计算a每列的平均值 C．a增加一行平均值 D．a增加一列平均值 3．在MATLAB命令行窗口输入下列命令： >> x=[1,2,3,4]; >> y=polyval(x,1); 则y的值为（）。 D A．5 B．8 C．24 D．10 4．设P是多项式系数向量，A为方阵，则函数polyval(P,A)与函数polyvalm(P,A)的值（）。D A．一个是标量，一个是方阵 B．都是标量 C．值相等 D．值不相等 5．在MATLAB命令行窗口输入下列命令： >> A=[1,0,-2]; >> x=roots(A); 则x(1)的值为（）。 C A．1 B．-2 C． D． 6．关于数据插值与曲线拟合，下列说法不正确的是（）。A A．3次样条方法的插值结果肯定比线性插值方法精度高。 B．插值函数是必须满足原始数据点坐标，而拟合函数则是整体最接近原始数据点，而不一定要必须经过原始数据点。 C．曲线拟合常常采用最小二乘原理，即要求拟合函数与原始数据的均方误差达到极小。 D．插值和拟合都是通过已知数据集来求取未知点的函数值。 二、填空题 1．设A=[1,2,3;10 20 30;4 5 6]，则sum(A)= ，median(A)= 。 [15 27 39]，[4 5 6[ 2．向量[2,0,-1]所代表的多项式是。2x2-1 3．为了求ax2+bx+c=0的根，相应的命令是（假定a、b、c已经赋值）。为了

实验一数据处理方法MATLAB实现

实验一数据处理方法的MATLAB实现 一、实验目的 学会在MATLAB环境下对已知的数据进行处理。 二、实验方法 1. 求取数据的最大值或最小值。 2. 求取向量的均值、标准方差和中间值。 3．在MATLAB环境下，对已知的数据分别进行曲线拟合和插值。 三、实验设备 1．586以上微机，16M以上内存，400M硬盘空间，2X CD-ROM 2．MATLAB5.3以上含CONTROL SYSTEM TOOLBOX。 四、实验内容 1．在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 2．在MATLAB环境下，选择合适的曲线拟合和插值方法，编写程序，对已知的数据分别进行曲线拟合和插值。 五、实验步骤 1. 在MATLAB环境下，将已知的数据存到数据文件mydat.mat中。 双击打开Matlab，在命令窗口（command window）中，输入一组数据：实验一数据处理方法的MATLAB实现 一、实验目的 学会在MATLAB环境下对已知的数据进行处理。 二、实验方法 1. 求取数据的最大值或最小值。 2. 求取向量的均值、标准方差和中间值。 3．在MATLAB环境下，对已知的数据分别进行曲线拟合和插值。 三、实验设备 1．586以上微机，16M以上内存，400M硬盘空间，2X CD-ROM 2．MATLAB5.3以上含CONTROL SYSTEM TOOLBOX。 四、实验内容

1．在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 2．在MATLAB环境下，选择合适的曲线拟合和插值方法，编写程序，对已知的数据分别进行曲线拟合和插值。 五、实验步骤 1. 在MATLAB环境下，将已知的数据存到数据文件mydat.mat中。 双击打开Matlab，在命令窗口（command window）中，输入一组数据： x=[1,4,2,81,23,45] x = 1 4 2 81 2 3 45 单击保存按钮，保存在Matlab指定目录（C:\Program Files\MATLAB71）下，文件名为“mydat.mat”。 2. 在MATLAB环境下，利用MATLAB控制系统工具箱中的函数直接求取数据的最大值或最小值，以及向量的均值、标准方差和中间值。 继续在命令窗口中输入命令： （1）求取最大值“max（a）”； >> max(x) ans = 81 （2）求取最小值“min（a）”； >> min(x) ans = 1 （3）求取均值“mean（a）”； >> mean(x) ans =

基于MATLAB的EXCEL数据计算与分析

基于MATLAB的EXCEL数据计算与分析 潜刘方 摘要：再怎么样希望先看摘要，阅读本文需要一定的MA TLAB基础知识，不需要excel相关知识。结合本人近期工作上的需要测量计算，想偷懒就选择了利用MATLAB偷懒，于是便有了本文。本文首先利用MA TLAB读取数据，计算，将数据写入excel，然后花了很大的精力来根据实际需要画图，最后将图保存在excel所在的文件夹下。这个m文件可谓花了我不少的时间和精力。最后根据m文件的不足（不能将图形输入到excel文档当中），进一步弥补这不足，就有了exlink（也叫excel link），在网上搜索了相关的知识，发现很多关于exlink 的培训，觉得实在可笑，所以就将exlink的使用写的比较详细，以供读者自行分析体会。关键字：MATLAB excel exlink 接口 一、前沿 MATLAB是一款应用在各个领域的数学软件，最初叫做矩阵实验室，专用于矩阵的运算，后来的版本再各个领域都得到了很好的应用，比如：通信、电力电子、电机控制、运动控制、计算机控制、自动控制，DSP数字信号处理。但是MATLAB对于数据的处理与可视化是很多软件所不能及的。 EXCEL作为办公必备软件，能对简单数据分析计算与作图分析，但是处理复杂数据显得力不从心，比如三维作图就无法利用EXCEL作出；EXCEL本身的函数远远没有MATLAB 多，MATLAB作为数据有其独特的优势，集成了很多数学函数，包括数据拟合差值等。MATLAB 可以从EXCEL中读取数据，经过相关运算之后又可以将数据写入EXCEL，假如需要重复性的对excel可以利用MATLAB编写函数，每次只要运行MATLAB程序就可以完成，大大节省时间和精力。 另外，MATLAB还有与EXCEL的接口，叫做EXLINK，运用这个接口可以在excel中完成MATLAB函数的调用，还能传送数据给MATLAB，从MATLAB当中读取数据，从MATLAB 当中读取图形，使用方便，操作简单。 二、基于MATLAB的数据分析 数据分析操作流程主要分为三步：第一步，从excel中读取数据；第二部：利用MATLAB 大量函数对数据分析处理；第三步：将分析结果写入excel中。在整个过程中，不需要打开excel软件，操作十分方便，每次操作唯一要做就是修改excel所在的目录及文件名。主要函数如下（具体使用方法可在MATLAB命令窗口输入help +函数名查看）：Xlsread 从excel中读数据 Xlswrite 向excel中邪数据 num2str 将数字转换为字符串 strncmp 字符串比较 polyfit 数据拟合 polyval 具体数值代入求值 plot 作图

第6章MATLAB数据分析与多项式计算_习题答案

第6章MATLAB数据分析与多项式计算 习题6 一、选择题 1．设A=[1,2,3,4,5;3,4,5,6,7]，则min(max(A))的值是（）。B A．1B．3C．5D．7 2．已知a为3×3矩阵，则运行mean(a)命令是（）。B A．计算a每行的平均值B．计算a每列的平均值 C．a增加一行平均值D．a增加一列平均值 3．在MATLAB命令行窗口输入下列命令： >>x=[1,2,3,4]; >>y=polyval(x,1); 则y的值为（）。D A．5B．8C．24D．10 4．设P是多项式系数向量，A为方阵，则函数polyval(P,A)与函数polyvalm(P,A)的值（）。D A．一个是标量，一个是方阵B．都是标量 C．值相等D．值不相等 5．在MATLAB命令行窗口输入下列命令： >>A=[1,0,-2]; >>x=roots(A); 则x(1)的值为（）。C A．1B．-2C．1.4142D．-1.4142 6．关于数据插值与曲线拟合，下列说法不正确的是（）。A A．3次样条方法的插值结果肯定比线性插值方法精度高。 B．插值函数是必须满足原始数据点坐标，而拟合函数则是整体最接近原始数据点，而不一定要必须经过原始数据点。 C．曲线拟合常常采用最小二乘原理，即要求拟合函数与原始数据的均方误差达到 极小。 D．插值和拟合都是通过已知数据集来求取未知点的函数值。 二、填空题 1．设A=[1,2,3;102030;456]，则sum(A)=，median(A)=。 [152739]，[456[ 2．向量[2,0,-1]所代表的多项式是。2x 2-1

2 3 ．为 了求a x +b x +c =0的根，相应的命令是（假定a 、b 、c 值）。为了 将求得的根代回方程进行验证，相应的命令是。 x=roots([a,b,c])，polyval([a,b,c],x) 4．如果被插值函数是一个单变量函插值，相应的MATLAB 函数 是。一维，interp1 5．求曲线拟合多项式系数的函数是，计算多项式在给定点上函数值的函数 是。polyfit ，polyval 三、应用题 1．利用M A T L A B 提供的r a n d n 函数生 成符合 正态分 布的10× 5随A ，进行如下 操作： （1）A 各列元素的均值和标。 （2）A 的最大元素和最小元素。 （3）求A 每行元素的和以及全部元素之和。 （4对A 的每列元素按升序、每行元素按降序排序。 第一题: (1): A=randn(10,5) B=mean(A) C=std(A) (2): mx=max(max(A)) mn=min(min(A)) (3): sm=sum(A,2) sz=sum(sum(A)) (4): [Y,I]=sort(A,1) [Z,J]=sort(A,2)； rot90(Z,1)'%旋转90度后，再转置便可得到每行按降序排列 2．已知多项式P1(x)=3x+2，P2(x)=5x 2-x+2，P3(x)=x 2-0.5，求： （1）P(x)=P 1(x)P2(x)P3(x)。 （2）P(x)=0的全部根。 （3）计算x i =0.2i(i=0，1，2，?，10)各点上的P(x i )。 第二题： (1): p1=[0,3,2]; p2=[5,-1,2]; 2