帮我设计一个IIR低通滤波器？

低通采样定理实验

11 实验目的

1．了解数字信号处理系统的一般构成；

2．掌握奈奎斯特抽样定理。

12 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

13 实验原理

一个典型的DSP系统除了数字信号处理部分外，还包括A/D和D/A两部分。这是因为自然界的信号，如声音、图像等大多是模拟信号，因此需要将其数字化后进行数字信号处理，模拟信号的数字化即称为A/D转换。数字信号处理后的数据可能需还原为模拟信号，这就需要进行D/A转换。一个仅包括A/D和D/A两部分的简化数字信号处理系统功能如图1所示。

A/D转换包括三个紧密相关的过程，即抽样、量化和编码。A/D转换中需解决的以下几个重要问题：抽样后输出信号中还有没有原始信号的信息？如果有能不能把它取出来？抽样频率应该如何选择？

奈奎斯特抽样定理（即低通信号的均匀抽样定理）告诉我们，一个频带限制在0至fx以内的低通信号x(t)，如果以fs≥2fx的抽样速率进行均匀抽样，则x(t)可以由抽样后的信号xs(t)完全地确定，即xs(t)包含有x(t)的成分，可以通过适当的低通滤波器不失真地恢复出x(t)。最小抽样速率fs=2fx称为奈奎斯特速率。

低通

译码

编码

量化

抽样

输入信号 样点输出 滤波输出

A/D（模数转换） D/A（数模转换）

图1 低通采样定理演示

为方便实现，实验中更换了一种表现形式，即抽样频率固定（10KHz），通过改变输入模拟信号的频率来展示低通抽样定理。我们可以通过研究抽样频率和模拟信号最高频率分量的频率之间的关系，来验证低通抽样定理。

14 实验内容

1．软件仿真实验：编写并调试MATLAB程序，分析有关参数，记录有关波形。

2．硬件实验：输入不同频率的正弦信号，观察采样时钟波形、输入信号波形、样点输出波形和滤波输出波形。

15 MATLAB参考程序和仿真内容

%%

%f—余弦信号的频率

% M—基2 FFT幂次数 N=2^M为采样点数，这样取值是为了便于作基2的FFT分析

%2 采样频率Fs

%%

function samples(f,Fs,M)

N=2^M; % fft点数=取样总点数

Ts=1/Fs; % 取样时间间隔

T=NTs; % 取样总时间=取样总点数取样时间间隔

n=0:N-1;

t=nTs;

Xn=cos(2fpit);

subplot(2,1,1);

stem(t,Xn);

axis([0 T 11min(Xn) 11max(Xn)]);

xlabel('t -->');

ylabel('Xn');

Xk=abs(fft(Xn,N));

subplot(2,1,2);

stem(n,Xk);

axis([0 N 11min(Xk) 11max(Xk)]);

xlabel('frequency -->');

ylabel('!Xk!');

%%

假如有一个1Hz的余弦信号y=cos(2πt)，对其用4Hz的采样频率进行采样，共采样32点，只需执行samples(1,4,5)，即可得到仿真结果。

软件仿真实验内容如下表所示：

仿真参数

f

Fs

Wo(计算)

Xn（图形）

Xk（图形）

（1，4，5）

另外记录图形，并标图号

（1，8，5）

（2，8，6）

自 选

16 硬件实验步骤

本实验箱采样频率fs固定为10KHz，低通滤波器的截止频率约为45KHz。

1、用低频信号源产生正弦信号，正弦信号源频率f自定，并将其接至2TP2（模拟输入）端，将示波器通道一探头接至2TP6（采样时钟）端观察采样时钟波形，示波器通道二探头接至2TP2观察并记录输入信号波形。

2、将示波器通道二探头接至2TP3观察并记录样点输出波形。

3、将示波器通道二探头接至2TP4观察并记录滤波输出波形。

4、根据采样定理，分f=fs /8、f=fs/4、f=fs/2等3种情况更改正弦信号频率，重复步骤2至步骤3。

5、用低频信号源产生方波信号，重复步骤1至步骤4。

17 思考题

1、 讨论在仿真实验中所计算的数字域频率Wo和Xk的图形中非零谱线位置之间的对应关系。

2、 讨论在仿真实验中自选参数的意义。

3、将在2TP2端加方波信号后的恢复波形，与相同频率的正弦信号的恢复波形相比，能够得出哪些结论？

2 FFT频谱分析实验

21 实验目的

1．通过实验加深对快速傅立叶变换（FFT）基本原理的理解。

2．了解FFT点数与频谱分辨率的关系，以及两种加长序列FFT与原序列FFT的关系。

22 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

23 实验原理

离散傅里叶变换（DFT）和卷积是信号处理中两个最基本也是最常用的运算，它们涉及到信号与系统的分析与综合这一广泛的信号处理领域。实际上卷积与DFT之间有着互通的联系：卷积可化为DFT来实现，其它的许多算法，如相关、滤波和谱估计等都可化为DFT来实现，DFT也可化为卷积来实现。

对N点序列x(n)，其DFT变换对定义为：

在DFT运算中包含大量的重复运算。FFT算法利用了蝶形因子WN的周期性和对称性，从而加快了运算的速度。FFT算法将长序列的DFT分解为短序列的DFT。N点的DFT先分解为2个N/2点的DFT，每个N/2点的DFT又分解为2个N/4点的DFT。按照此规律，最小变换的点数即所谓的“基数（radix）。”因此，基数为2的FFT算法的最小变换（或称蝶形）是2点DFT。一般地，对N点FFT，对应于N个输入样值，有N个频域样值与之对应。一般而言，FFT算法可以分为时间抽取（DIT）FFT和频率抽取（DIF）两大类。

在实际计算中，可以采用在原来序列后面补0的加长方法来提高FFT的分辨率；可以采用在原来序列后面重复的加长方法来增加FFT的幅度。

24 实验内容

1．软件仿真实验：分别观察并记录正弦序列、方波序列及改变FFT的点数后的频谱；分别观察并记录正弦序列、方波序列及2种加长序列等信号的频谱。

2．硬件实验：分别观察并记录正弦信号、方波信号及改变FFT的点数后的频谱。

25 MATLAB参考程序和仿真内容

%%

function[x]=ffts(mode,M)

Nfft=2^M;

x=zeros(1,Nfft); %定义一个长度为Nfft的一维全0数组

if mode= =1 for n=0:Nfft-1 x(n+1)=sin(2pin/Nfft); end

end %定义一个长度为Nfft的单周期正弦序列

if mode= =2 for n=0:Nfft-1 x(n+1)=sin(4pin/Nfft); end

end %定义一个长度为Nfft的双周期正弦序列

if mode= =3 for n=0:Nfft/2-1 x(n+1)=sin(4pin/Nfft); end

end %定义一个长度为Nfft/2的正弦序列，后面一半为0序列。

if mode= =4 for n=0:Nfft-1 x(n+1)=square(2pin/Nfft); end

end

if mode= =5 for n=0:Nfft-1 x(n+1)=square(2pin/Nfft); end

end

if mode= =6 for n=0:Nfft/2-1 x(n+1)=square(4pin/Nfft); end

end

n=0:Nfft-1;

subplot(2,1,1);

stem(n,x);

axis([0 Nfft-1 11min(x) 11max(x)]);

xlabel('Points-->');

ylabel('x(n)');

X=abs(fft(x,Nfft));

subplot(2,1,2);

stem(n,X);

axis([0 Nfft-1 11min(X) 11max(X)]);

xlabel('frequency-->');

ylabel('!X(k)!');

%%

假设需观察方波信号的频谱，对一个周期的方波信号作32点的FFT，则只需在MATLAB的命令窗口下键入：[x]=ffts(21,5) ，程序进行模拟，并且输出FFT的结果。

关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个以上的例子进行记录。例如要观察后面补0的加长方法来提高FFT的分辨率的现象，可以仿真ffts(4,5)和ffts(6,6)两个例子。

26 硬件实验步骤

1．将低频信号源输出加到实验箱模拟通道1输入端，将示波器探头接至模拟通道1输出端。

2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“FFT频谱分析”子菜单，出现显示FFT频谱分析功能提示信息的窗口。

3．用低频信号产生器产生一个1KHz的正弦信号。

4．选择FFT频谱分析与显示的点数为64点，开始进行FFT运算。此后，计算机将周期性地取回DSP运算后的FFT数据并绘图显示

5．改信号源频率，观察并记录频谱图的变化。

6．选择FFT的点数为128点，观察并记录频谱图的变化。

7．更改正弦信号的频率，重复步骤4 ~步骤6。

8．用低频信号产生器产生一个1KHz的方波信号，重复步骤4 ~步骤7。注意：应根据实验箱采样频率fs为10KHz和方波信号的频带宽度选择方波信号的频率。

本硬件实验要进行两种信号，每个信号两种频率，每个信号两种点数等共8次具体实验内容，性质能够体现实验要求的4个以上的例子进行记录。

27 思考题

1．对同一个信号，不同点数FFT观察到的频谱图有何区别？

2．序列加长后FFT与原序列FFT的关系是什么，试推导其中一种关系。

3．用傅立叶级数理论，试说明正弦信号频谱和方波信号频谱之间的关系。

3 IIR滤波器设计实验

31 实验目的

1．通过实验加深对IIR滤波器基本原理的理解。

2．学习编写IIR滤波器的MATLAB仿真程序。

32 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

33 实验原理

IIR滤波器有以下几个特点：

1．IIR数字滤波器的系统函数可以写成封闭函数的形式。

2．IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。

3．IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

4．IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。

在MATLAB下设计IIR滤波器可使用Butterworth函数设计出巴特沃斯滤波器，使用Cheby1函数设计出契比雪夫I型滤波器，使用Cheby2设计出契比雪夫II型滤波器，使用ellipord函数设计出椭圆滤波器。下面主要介绍前两个函数的使用。

与FIR滤波器的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。

一、巴特沃斯IIR滤波器的设计

在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。

Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。

butter函数的用法为：

[b,a]=butter(n,Wn,/ftype/)

其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为：

[n,Wn]= buttord(Wp,Ws,Rp,Rs)

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则：

1．高通滤波器：Wp和Ws为一元矢量且Wp>Ws；

2．低通滤波器：Wp和Ws为一元矢量且Wp<Ws；

3．带通滤波器：Wp和Ws为二元矢量且Wp<Ws，如Wp=[02,07],Ws=[01,08];

4．带阻滤波器：Wp和Ws为二元矢量且Wp>Ws，如Wp=[01,08],Ws=[02,07]。

二、契比雪夫I型IIR滤波器的设计

在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。

cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。

cheby1函数的用法为：

[b,a]=cheby1(n,Rp,Wn,/ftype/)

在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。

cheblord函数的用法为：

[n,Wn]=cheblord(Wp,Ws,Rp,Rs)

其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。

34 实验内容

1．软件仿真实验：编写并调试MATLAB程序，选择不同形式，不同类型的4种滤波器进行仿真，记录幅频和相频特性，对比巴特沃斯滤波器和契比雪夫滤波器。

2．硬件实验：设计IIR滤波器，在计算机上观察冲激响应、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。

35 MATLAB参考程序和仿真内容

%%

%mode: 1--巴特沃斯低通；2--巴特沃斯高通；3--巴特沃斯带通；4--巴特沃斯带阻

% 5--契比雪夫低通；6--契比雪夫高通；7--契比雪夫带通；8--契比雪夫带阻

%fp1,fp2： 通带截止频率，当高通或低通时只有fp1有效

%fs1, fs2： 阻带截止频率，当高通或低通时只有fs1有效

%rp: 通带波纹系数

%as: 阻带衰减系数

%sample: 采样率

%h: 返回设计好的滤波器系数

%%

function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)

wp1=2fp1/sample;wp2=2fp2/sample;

ws1=2fs1/sample;ws2=2fs2/sample;

%得到巴特沃斯滤波器的最小阶数N和3bd频率wn

if mode<3[N,wn]=buttord(wp1,ws1,rp,as);

elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as);

%得到契比雪夫滤波器的最小阶数N和3bd频率wn

elseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as);

else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as);

end

%得到滤波器系数的分子b和分母a

if mode= =1[b,a]=butter(N,wn);end

if mode= =2[b,a]=butter(N,wn,/high/);end

if mode= =3[b,a]=butter(N,wn);end

if mode= =4[b,a]=butter(N,wn,/stop/);end

if mode= =5[b,a]=cheby1(N,rp,wn);end

if mode= =6[b,a]=cheby1(N,rp,wn,/high/);end

if mode= =7[b,a]=cheby1(N,rp,wn);end

if mode= =8[b,a]=cheby1(N,rp,wn,/stop/);end

set(gcf,/menubar/,menubar);

freq_response=freqz(b,a);

magnitude=20log10(abs(freq_response));

m=0:511;

f=msample/(2511);

subplot(3,1,1);plot(f,magnitude);grid; %幅频特性

axis([0 sample/2 11min(magnitude) 11max(magnitude)]);

ylabel('Magnitude');xlabel('Frequency-->');

phase=angle(freq_response);

subplot(3,1,2);plot(f,phase);grid; %相频特性

axis([0 sample/2 11min(phase) 11max(phase)]);

ylabel('Phase');xlabel('Frequency-->');

h=impz(b,a,32); %32点的单位函数响应

t=1:32;

subplot(3,1,3);stem(t,h);grid;

axis([0 32 12min(h) 11max(h)]);

ylabel('h(n)');xlabel('n-->');

%%

假设需设计一个巴特沃斯低通IIR滤波器，通带截止频率为2KHz，阻带截止频率为3KHz，通带波纹系数为1，阻带衰减系数为20，采样频率为10KHz，则只需在MATLAB的命令窗口下键入：

[b,a]=iirfilt(1,2000,3000,2400,2600,1,20,10000)

程序进行模拟，并且按照如下顺序输出数字滤波器系统函数

的系数

b= b0 b1 ……bn

a= a0 a1 ……an

关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个例子进行记录，系统函数只要记录系统的阶数。

36 硬件实验步骤

1．根据实验箱采样频率fs为10KHz的条件，用低频信号发生器产生一个频率合适的低频正弦信号，将其加到实验箱模拟通道1输入端，将示波器通道1探头接至模拟通道1输入端，通道2探头接至模拟通道2输出端。

2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“IIR滤波器”子菜单，出现提示信息。

3．输入滤波器类型、滤波器截止频率等参数后，分别点击“幅频特性”和“相频特性”按钮，在窗口右侧观察IIR滤波器的幅频特性和相频特性。此时提示信息将消失，如需查看提示信息，可点击“设计说明”按钮。

4．点击“下载实现”按钮，IIR滤波器开始工作，此时窗口右侧将显示IIR滤波器的幅频特性。

5．根据输入滤波器类型，更改低频信号源的频率，观察示波器上输入输出波形幅度的变化情况，测量IIR滤波器的幅频响应特性，看其是否与设计的幅频特性一致。

6．更改滤波器类型、滤波器截止频率等参数（共4种），重复步骤3至步骤5。所选择的例子参数最好和MATLAB仿真程序的例子一样。

7．用低频信号产生器产生一个500Hz的方波信号，分别设计3种滤波器，完成如下表要求的功能，并且记录参数和波形。

功 能

滤波器类型

参 数

输出波形

fp1

fp2

fs1

fs2

通过3次及以下次数的谐波

另外记录图形，并标图号

滤除5次及以下次数的谐波

通过3次到5次的谐波

37 思考题

1．在实验箱采样频率fs固定为10KHz的条件下，要观察方波信号频带宽度内的各个谐波分量，方波信号的频率最高不能超过多少，为什么？

2．硬件实验内容7中输出信号各个谐波分量，与原来方波信号同样谐波分量相比，有没有发生失真？主要发生了什么类型的失真？为什么？

4 窗函数法FIR滤波器设计实验

41 实验目的

1．通过实验加深对FIR滤波器基本原理的理解。

2．学习使用窗函数法设计FIR滤波器，了解窗函数的形式和长度对滤波器性能的影响。

42 实验仪器

1．YBLD智能综合信号源测试仪 1台

2．双踪示波器 1台

3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台

4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台

43 实验原理

数字滤波器的设计是数字信号处理中的一个重要内容。数字滤波器设计包括FIR（有限单位脉冲响应）滤波器与IIR（无限单位脉冲响应）滤波器两种。

与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数H（z）为：

H（z）是z－1的N－1次多项式，它在z平面上有N－1个零点，原点z=0是N－1阶重极点，因此H（z）是永远稳定的。稳定和线性相位特性是FIR滤波器突出的优点。

FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H( )满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本实验介绍窗函数法的FIR滤波器设计。

窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标准响应的高通、低通、带通和带阻FIR滤波器。

一、firl函数的使用

在MATLAB下设计标准响应FIR滤波器可使用firl函数。firl函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为：

b=firl(n,Wn,/ftype/,Window)

各个参数的含义如下：

b—滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：b(z)=b(1)+b(2)z－1+…+b(n+1)z－n。

n—滤波器阶数。

Wn—截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时，Wn=[W1 W2],W1≤ω≤W2。

ftype—当指定ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。

Window—窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数n+1。

二、窗函数的使用

在MATLAB下，这些窗函数分别为：

1．矩形窗：w=boxcar(n)，产生一个n点的矩形窗函数。

2．三角窗：w=triang(n)，产生一个n点的三角窗函数。

当n为奇数时，三角窗系数为w(k)=

当n为偶数时，三角窗系数为w(k)=

3．巴特利特窗：w=Bartlett(n)，产生一个n点的巴特利特窗函数。

巴特利特窗系数为w(k)=

巴特利特窗与三角窗非常相似。巴特利特窗在取样点1和n上总以零结束，而三角窗在这些点上并不为零。实际上，当n为奇数时bartlett(n)的中心n－2个点等效于triang(n－2)。

4．汉明窗：w=hamming(n)，产生一个n点的汉明窗函数。

汉明窗系数为w(k+1)=054－046cos（ ） k=0,…,n－1

5．汉宁窗：w=hanning(n)，产生一个n点的汉宁窗函数。

汉宁窗系数为w(k)=05[1－cos( )] k=1,…,n

6．布莱克曼窗：w=Blackman(n),产生一个n点的布莱克曼窗函数。

布莱克曼窗系数为w(k)=042－05cos(2π )+08cos（4π ）] k=1,…,n

与等长度的汉明窗和汉宁窗相比，布莱克曼窗的主瓣稍宽，旁瓣稍低。

7．凯泽窗：w=Kaiser(n,beta)，产生一个n点的凯泽窗数，其中beta为影响窗函数旁瓣的β参数，其最小的旁瓣抑制α与β的关系为:

01102(α－087) α>50

β= 05842(α－21)04+007886(α－21) 21≤α≤50

0 α<21

增加β可使主瓣变宽，旁瓣的幅度降低。

8．契比雪夫窗：w=chebwin(n,r)产生一个n点的契比雪夫窗函数。其傅里叶变换后的旁瓣波纹低于主瓣r个db数。

44 实验内容

1．软件仿真实验：编写并调试MATLAB程序，观察不同窗，不同类型滤波器不同点数等共4种FIR滤波器的h(n)，并记录幅频特性和相频特性。

2．硬件实验：用窗函数法设计标准响应的FIR滤波器，在计算机上观察窗函数幅频特性、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。

45 MATLAB参考程序和仿真内容

%%

%mode: 模式（1--高通；2--低通；3--带通；4--带阻）

%n: 阶数，加窗的点数为阶数加1

%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率

%fs: 带通和带阻时指示上限频率

%window:加窗（1--矩形窗；2--三角窗；3--巴特利特窗；4--汉明窗；

% 5--汉宁窗；6--布莱克曼窗；7--凯泽窗；8--契比雪夫窗）

%r: 代表加chebyshev窗的r值和加kaiser窗时的beta值

%sample: 采样率

%h: 返回设计好的FIR滤波器系数

%%

%mode: 模式（1--高通；2--低通；3--带通；4--带阻）

%n: 阶数，加窗的点数为阶数加1

%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率

%fs:

题目：利用DSP的FIR滤波器设计

数字处理器（DSP）有很强的数据处理能力，它在高速数字信号处理领域有广泛的应用，例如数字滤波、音频处理、图像处理等。相对于模拟滤波器，数字滤波器没有漂移，能够处理低频信号，频率响应特性可做成非常接近于理想的特性，且精度可以达到很高，容易集成等。使用可编程的DSP芯片实现数字滤波可以通过修改滤波器的参数十分方便地改变滤波器的特性，下面主要说明利用TMS320VC54x DSP芯片设计实现FIR数字滤波器。

设计目的意义

一个实际的应用系统中,总存在各种干扰，所以在系统设计中,滤波器的好坏将直接影响系统的性能。使用DSP进行数字处理，可以对一个具有噪声和信号的混合信号源进行采样,再经过数字滤波,滤除噪声,就可以提取有用信号了。所以说，数字滤波器是DSP最基本的应用领域，熟悉基于DSP的数字滤波器能为DSP应用系统开发提供良好的基础。

技术指标

1、数字滤波器的频率参数主要有：①通带截频：为通带与过渡带的边界点，在该点信号增益下降到规定的下限。②阻带截频：为阻带与过渡带的边界点，在该点信号衰耗下降到规定的下限。③转折频率：为信号功率衰减到1/2(约3dB)时的频率，在很多情况下，也常以fc作为通带或阻带截频。④当电路没有损耗时，固有频率：就是其谐振频率，复杂电路往往有多个固有频率。

2、增益与衰耗

滤波器在通带内的增益并非常数。①对低通滤波器通带增益，一般指ω=0时的增益；高通指ω→∞时的增益；带通则指中心频率处的增益。②对带阻滤波器，应给出阻带衰耗，衰耗定义为增益的倒数。③通带增益变化量指通带内各点增益的最大变化量，如果通带增益变化量以dB为单位，则指增益dB值的变化量。

3、阻尼系数与品质因数

阻尼系数α是表征滤波器对角频率为ω0信号的阻尼作用，是滤波器中表示能量衰耗的一项指标，它是与传递函数的极点实部大小相关的一项系数。

4、灵敏度

滤波电路由许多元件构成，每个元件参数值的变化都会影响滤波器的性能。

5、群时延函数

在滤波器设计中，常用群时延函数评价信号经滤波后相位失真程度。

以上的几个技术指标是一般滤波器的特性，但在实际应用中，数字滤波器通常用来实现选频 *** 作，因此在利用DSP实现数字滤波器设计中要求的技术指标主要为在频域中给出的幅频响应和相频响应。如下图所示

幅频响应和相频响应特性曲线

对于幅频响应，它的含义是信号通过系统之后的输出信号的幅度与它输入时的信号的幅度的比值，一般以分贝值表示。对于相频响应，含义是信号通过系统之后的输出信号的相位与它输入时的信号的相位之差，在运用线性相频响应指标进行滤波器设计具有如下优点：①只包含实数算法，不涉及复数运算；②不存在延迟失真，只有固定数量的延迟；③可以采用FFT算法，从而提高运行效率；④由于FIR滤波器的单位脉冲响应是有限长序列，故FIR滤波器没有不稳定的问题,且误差较小。

基本原理

利用DSP实现FIR滤波器的设计方法主要有窗函数法和频率抽样法，其中窗函数法是基本的设计方法，这里采用窗函数法设计FIR滤波器。设希望得到的滤波器理想响应为 ,那么FIR滤波器的设计就在于寻找一个传递函数

去逼进 ，设

这里 就是傅立叶级数的系数。在这种逼近中，最直接的一种方法就是从单位脉冲响应 入手，使 逼近理想的单位脉冲响应 。由于 是一个无限长序列，因此，最简单的方法就是对 做截尾处理，即得到一个近似的传递函数

上式中，Q就是最终确定FIR滤波器的阶数，Q越大，近似程度就越高。对 截尾，实际上就是对 乘上一个矩形窗口 ,即

令z= ,则

其脉冲响应系数为 ， ，…， ， ， ，…， ， 。为使 具有因果性，延时Q个样值，可得：

令n+Q=k，上式成为

令 ,N=2Q,得

式中， 是脉冲响应系数，这里 …, ,…, 。

一般来说,FIR数字滤波器输出 的Z变换形式 与输入 的Z变换形式之间的关系如下:

实现结构如下图所示：

Z变换结构图

从上面的Z变换和结构图可以很容易得出FIR滤波器的差分方程表示形式，即对上式进行反Z变换得：

上式为FIR数字滤波器的时域表示方法,其中x(n)是在时间n的滤波器的输入抽样值，根据上式即可对滤波器进行设计。

硬件设计

1、DSP芯片

根据设计原理，实现的核心器件采用美国德州仪器公司生产的低功耗定点数字信号处理器芯片TMS320C5402。选择该芯片主要是因为它是目前最常用的低成本DSP芯片，而且包括以下主要特点：

⑴运算速度快，最快可达532MIPS；

⑵多总线结构，片内共有8 条总线（1条程序存储器总线、3条数据存储总线和4条地址总线）；

⑶CPU采用冯 诺依曼并行结构设计，使其能在一条指令周期内，高速地完成多项算术运算；

⑷片内集成了4K×16bitROM和16K×16bit的双存取RAM；

⑸丰富的片上外围电路（通用I/O 引脚，定时器，时钟发生器， HPI 接口，多通道缓冲串行口McBSP）使其与外部接口方便；

⑹33V I/O电压，18V核点压，工作电流平均值为75mA，其中核45mA，I/O约30mA；

⑺144脚BGA封装，使体积减少，功耗降低。

2、AD和DA电路

在本数字滤波器系统中选择了TI公司的TLV1570芯片作为模数转换器件，8通道10位27到55 V低电压模数转换芯片。TLVl570在3V电压下的采样频率为625KSPS，输入信号最高频率不能超过300K。

由于模数转换选择了10位器件，为了简化程序代码，减少DSP 的运算工作量，在本数字滤波器系统中选择了TI公司的TLV5608芯片，它是一款8通道10位27到55V低电压数模转换芯片。

3、电源电路

根据DSP芯片工作的电压电流需求，及芯片采用双电源供电对加电顺序的要求，考虑使用TI公司的电源转换芯片TPS73HD318,其输出电压为一路33V、一路18V，每路电源的最大输出电流为750mA，能满足本系统的供电需求。而且TPS73xx具有非常低的静态电流，能使稳压器输出稳定。

4、时钟电路

C54xx系列的时钟端子为X1和X2/CLKIN，采用无源晶振提供时钟信号，由于DSP有一组端子可以用来调整其工作频率的高低，故对晶振频率大小的选定没有特别的要求，这里选用10Mhz的晶振。

5、复位电路

为了克服DSP系统因时钟频率较高导致在运行时可能发生的干扰和被干扰的现象，最好是使用具有监视（Watchdog）功能的自动复位电路,于是采用专门的自动复位芯片MAX706。MAX706的电源为31V～50V,低电平复位输出，复位门限为308V。

6、未用端子处理

根据使用DSP芯片的相关原则，以及芯片手册具体决定未用端子是接上拉电阻还是悬空。

7、基于上述的各部分电路组成，可以得出DSP数字滤波器的整体硬件电路连线图，如下所示

程序设计

1、设计思路

⑴在DSP进行数字滤波运算前首先要进行初始化，只有正确设置了DSP的初始状态才能保证芯片能正常运行。本系统主要进行以下两方面的初始化：

①寄存器初始化：状态寄存器ST0、状态寄存器ST1、处理器模式控制寄存器PMST、软件等待状态寄存器SWWSR、组交换控制寄存器BSCR和时钟模式寄存器等。

②中断矢量表初始化：根据DSP芯片对各中断矢量的设置位置编写一个子程序；设置PMST控制寄存器；连接时将矢量表重定位到IPTR指定的地址。

⑵其次就是FIR 数字滤波的子程序设计，主要步骤如下：

①查询SPCR11寄存器的第二位，当为1时说明read ready，将DRR11的值读入AR3所指向的地址，该值为最新的采样值。

②将最新的采样值减去200h，然后AR3的值减1。

③执行MAC指令。

④将累加器的值送给变量Y，并将Y加上200h。

⑤查询SPCR20寄存器的第二位，当为1时说明writeready，将Y值赋给DXR10，该值为滤波器输出值。

⑥循环执行上面步骤。

2、程序流程图

依据上述程序设计思路可以得到利用DSP实现FIR滤波器设计的程序流程图，如下

3、程序代码

由于初始化程序部分过于庞大繁杂，这里只给出用MAC指令编程实现FIR低通滤波器的程序片断：

FILT_task1

LD Store_SICX,A

STLM A,ar4

STM #1,ar0 ；间址

STM #28,bk

LD DEM_Out,A

STL A,ar4+% ；输入信号:实部

STM #Coef_Tab1,ar5 ；滤波器实部系数地址

LD #0,A

STM #27,brc

RPTB SICXU-1

MAC AR4+0%,AR5+,A

SICXU LD A,-16,A ；低通滤波结果

LD C7FFF,B

MIN A

NEG B

MAX A

STL A,DEM_Out

LDM AR4,A

STL A,Store_SICX

RET

Coef_Tab1

word 100 ；h(0)

word 7 ；h(1)

… ；脉冲响应系数

word -248

word -71 ；h(N-1)

end

总结

通过利用DSP的FIR滤波器设计，对DSP芯片的使用，以及利用DSP芯片组成的基本系统的相关电路有了比较深的认识。熟悉DSP芯片的系统设计和应用开发流程，并利用图书馆、网络、询问同学等方式查找资料和解决相关的难题，这是最基础的工作，也是最关键的步骤。这样做可以培养自己的动手解决问题的能力和独立思考的处事方法，使自己具有技术人员的气质和工作态度，为将来的就业增加优势。

数字滤波器是DSP的典型应用，学会了有助于触类旁通，利于进一步的学习研究，能做到理解其他基于DSP的系统的功能和工作原理。掌握了基于DSP的应用开发，开阔了视野，增长了知识，是进入现代数字信号处理领域重要技能，乃至大规模集成电路的开发也是会用到的基础，今后要予以重视并积极努力去学习。

MATLAB设计模拟带通滤波器

参数自己改一下就可以了

cheb1

% wp1=045pi;wp2=065pi;ws1=03pi;ws2=075pi;Rp=1;Rs=40

% =============双线型变换法=========================================

wp1=045pi; wp2=065pi;

ws1=03pi; ws2=075pi;

Rp=1; Rs=40;

Wp1=tan(wp1/2); Wp2=tan(wp2/2);

Ws1=tan(ws1/2); Ws2=tan(ws2/2);

BW=Wp2-Wp1; W0=Wp1Wp2; W00=sqrt(W0);

WP=1; WS=WP(W0^2-Ws1^2)/(Ws1BW);

[N,Wn]=cheb1ord(WP,WS,Rp,Rs,'s');

[B,A]=cheby1(N,Rp,Wn,'s');

[BT,AT]=lp2bp(B,A,W00,BW);

[num,den]=bilinear(BT,AT,05);

[h,omega]=freqz(num,den,64);

subplot(2,2,1);stem(omega/pi,abs(h));

xlabel('\omega/\pi');ylabel('|H(z)|');

subplot(2,2,2);stem(omega/pi,20log10(abs(h)));

xlabel('\omega/\pi');ylabel('增益dB');

% =============直接法=================================

wp1=045pi; wp2=065pi;

ws1=03pi; ws2=075pi;

Rp=1; Rs=40;

Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];

[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs);

[B,A]=cheby1(N,Rp,Wn);

[h,omega]=freqz(B,A,64);

subplot(2,2,3);stem(omega/pi,abs(h));

xlabel('\omega/\pi');ylabel('|H(z)|');

subplot(2,2,4);stem(omega/pi,20log10(abs(h)));

xlabel('\omega/\pi');ylabel('增益dB');

%cheby2%

% wp1=045pi;wp2=065pi;ws1=03pi;ws2=075pi;Rp=1;Rs=40

% =============双线型变换法=========================================

wp1=045pi; wp2=065pi;

ws1=03pi; ws2=075pi;

Rp=1; Rs=40;

Wp1=tan(wp1/2); Wp2=tan(wp2/2);

Ws1=tan(ws1/2); Ws2=tan(ws2/2);

BW=Wp2-Wp1; W0=Wp1Wp2; W00=sqrt(W0);

WP=1; WS=WP(W0^2-Ws1^2)/(Ws1BW);

[N,Wn]=cheb2ord(WP,WS,Rp,Rs,'s');

[B,A]=cheby2(N,Rs,Wn,'s');

[BT,AT]=lp2bp(B,A,W00,BW);

[num,den]=bilinear(BT,AT,05);

[h,omega]=freqz(num,den,64);

subplot(2,2,1);stem(omega/pi,abs(h));

xlabel('\omega/\pi');ylabel('|H(z)|');

subplot(2,2,2);stem(omega/pi,20log10(abs(h)));

axis([0 1 -100 0]);xlabel('\omega/\pi');ylabel('增益dB');

% =============直接法=================================

wp1=045pi; wp2=065pi;

ws1=03pi; ws2=075pi;

Rp=1; Rs=40;

Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];

[N,Wn]=cheb2ord(Wp,Ws,Rp,Rs);

[B,A]=cheby2(N,Rs,Wn);

[h,omega]=freqz(B,A,64);

subplot(2,2,3);stem(omega/pi,abs(h));

xlabel('\omega/\pi');ylabel('|H(z)|');

subplot(2,2,4);stem(omega/pi,20log10(abs(h)));

axis([0 1 -100 0]);xlabel('\omega/\pi');ylabel('增益dB');

%椭圆%

% wp1=045pi;wp2=065pi;ws1=03pi;ws2=075pi;Rp=1;Rs=40

% =============双线型变换法=========================================

wp1=045pi; wp2=065pi;

ws1=03pi; ws2=075pi;

Rp=1; Rs=40;

Wp1=tan(wp1/2); Wp2=tan(wp2/2);

Ws1=tan(ws1/2); Ws2=tan(ws2/2);

BW=Wp2-Wp1; W0=Wp1Wp2; W00=sqrt(W0);

WP=1; WS=WP(W0^2-Ws1^2)/(Ws1BW);

[N,Wn]=ellipord(WP,WS,Rp,Rs,'s');

[B,A]=ellip(N,Rp,Rs,Wn,'s');

[BT,AT]=lp2bp(B,A,W00,BW);

[num,den]=bilinear(BT,AT,05);

[h,omega]=freqz(num,den,64);

subplot(2,2,1);stem(omega/pi,abs(h));grid;

xlabel('\omega/\pi');ylabel('|H(z)|');

subplot(2,2,2);stem(omega/pi,20log10(abs(h)));grid;

xlabel('\omega/\pi');ylabel('增益dB');

% =============直接法=================================

wp1=045pi; wp2=065pi;

ws1=03pi; ws2=075pi;

Rp=1; Rs=40;

Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];

[N,Wn]=ellipord(Wp,Ws,Rp,Rs);

[B,A]=ellip(N,Rp,Rs,Wn);

[h,omega]=freqz(B,A,64);

subplot(2,2,3);stem(omega/pi,abs(h));grid;

xlabel('\omega/\pi');ylabel('|H(z)|');

subplot(2,2,4);stem(omega/pi,20log10(abs(h)));grid;

xlabel('\omega/\pi');ylabel('增益dB');

回答者： dear漫儿 - 初学弟子 一级 5-6 19:37

上网搜了一下，还没发现把mean shift原理说的很清楚的文章，我就来写一段了：

首先你要知道这是PAMI 2003的一篇文章，非常经典的哦，Mean Shift: A Robust Approach Toward Feature Space Analysis。当然原文有17页，都是一些复杂的公式。。。。。

Mean shift主要用在图像平滑和图像分割（那个跟踪我现在还不清楚），先介绍一下平滑的原理：

输入是一个5维的空间，2维的（x，y）地理坐标，3维的（L，u，v）的颜色空间坐标，当然你原理也可以改写成rgb色彩空间或者是纹理特征空间。

先介绍一下核函数，有uniform的，也有高斯的核函数，不管是哪个的，其基本思想如下：简单的平滑算法用一个模板平均一下，对所有的像素，利用周围的像素平均一下就完事了，

这个mean shift的是基于概率密度分布来的，而且是一种无参的取样。有参的取样就是假设所有的样本服从一个有参数的概率分布函数，比如说泊松分布，正态分布等等，高中生都知道概率公式里面是有参数的，在说一下特征空间是一个5维的空间，距离用欧几里德空间就可以了，至少代码里就是这样实现的，而本文的无参取样是这样的：在特征空间里有3维的窗口（想象一下2维空间的窗口），对于一个特征空间的点，对应一个5维的向量，可以计算该点的一个密度函数，如果是有参的直接带入该点的坐标就可以求出概率密度了，基于窗函数的思想就是考虑它邻近窗口里的点对它的贡献， 它假设密度会往密集一点的地方转移，算出移动之后的一个5维坐标，该坐标并会稳定，迭代了几次之后，稳定的地方是modes。这样每一个像素点都对应一个这么一个modes，用该点的后3维的值就是平滑的结果了，当然在算每个点的时候，有些地方可能重复计算了，有兴趣的化你可以参考一下源代码，确实是可以优化的。总结一下mean shift的平滑原理就是在特征空间中向密度更高的地方shift（转移）。

其次是怎么利用mean shift分割图像

先对图像进行平滑，

第2步利用平滑结果建立区域邻接矩阵或者区域邻接链表，就是在特征空间比较近的二间在2维的图像平面也比较接近的像素算成一个区域，这样就对应一个区域的邻接链表，记录每个像素点的label值。当然代码中有一个传递凸胞的计算，合并2个表面张力很接近的相邻区域，这个我还没想怎么明白，希望比较清楚的朋友讲一讲。最后还有一个合并面积较小的区域的 *** 作，一个区域不是对应一个modes值嘛，在待合并的较小的那个区域中，寻找所有的邻接区域，找到距离最小的那个区域，合并到那个区域就ok了。

void msImageProcessor::Filter(int sigmaS, float sigmaR, SpeedUpLevel speedUpLevel)这个是平滑 *** 作，sigmaS是2维的平面窗口的窗函数，sigmaR是颜色空间的窗函数，最后一个参数表示是否加速算法。

void msImageProcessor::FuseRegions(float sigmaS, int minRegion)这个就是合并较小区域的一个函数。

void msImageProcessor::Segment(int sigmaS, float sigmaR, int minRegion, SpeedUpLevel speedUpLevel)这个是分割函数，里面包括了平滑功能。再详细的我也不说了。。。。。

最后是怎么调用问题：void CImageProcessing::mydebug()

{

int x,y;

int width_ = imageWidth;

int height_ = imageHeight;

unsigned char tmp = new unsigned char[width_ height_ 3];

unsigned char dst = tmp;

for (x = 0;x < imageHeight; x ++)

for(y = 0;y < imageWidth; y++)

{

//Color tmp = (imageData_RGB)(y, x);

//uchar data = &CV_IMAGE_ELEM(ip,uchar,x,y3) ;

// (dst++) = data[2] ;

// (dst++) = data[1] ;

// (dst++) = data[0] ;

(dst++) = (imageData_RGB)(y, x)channel[0];

(dst++) = (imageData_RGB)(y, x)channel[1];

(dst++) = (imageData_RGB)(y, x)channel[2];

}

cbgImage_->SetImageFromRGB( tmp, width_, height_, true);

delete []tmp;

msImageProcessor iProc = new msImageProcessor();

iProc->DefineImage(cbgImage_->im_, COLOR, height_, width_);

iProc->SetSpeedThreshold(01);

iProc->Segment(5,8,10,NO_SPEEDUP);

// iProc->Filter(5,8,NO_SPEEDUP);

// iProc->FuseRegions(mWinSize,mAreaSize);//°ë¾¶ºÍ×îСÇøÓòÃæ»ý

} 至于怎么调试编译我看我是说不清楚了。。。。。。。。。。。。。。。。。。。。。。

推荐参考：http://hibaiducom/zzf378139208/blog/item/33a47c35b7ed95b9d1a2d304html

http://hibaiducom/zzf378139208/blog/item/ee395c115019e80f203f2e85html

SciPy提供了firwin用窗函数设计低通滤波器，firwin的调用形式如下：

firwin(N, cutoff, width=None, window='hamming')

其中N为滤波器的长度；cutoff为以正规化的频率；window为所使用的窗函数。

冯强强1，2　温明明1，2　吴衡1，2

（1广州海洋地质调查局 广州 510760；2国土资源部海底矿产资源重点实验室 广州 510760）

第一作者简介：冯强强（1987—），男，硕士，助理工程师，现在从事海上地球物理探测方法的研究工作。

摘要 本文选取了两条不同的单道地震接收电缆（AAE接收电缆、GEO接收电缆）在野外采集到的原始数据，并利用快速傅立叶变换（FFT）以及S变换等方法，从数据本身、频率域、时频域方面进行了对比分析。结果表明接收段较短、且内部检波器组合更紧致的AAE电缆具有较高的浅层分辨率，更适合于浅水作业，而接收段较长，检波器较多，检波器间距稍大的GEO电缆可以采集到高信噪比的资料，更适合于深水作业。

关键词 单道地震 S 变换 时频域 检波器

1 引言

单道地震是解决海上工程地质问题以及海洋区域地质调查的一种重要手段［1，2］，相对于使用全波形反射的多道地震调查，单道地震对震源和检波器的偏移距离相对于水深和对反射倾角的要求很小，另外，由于其施工方法简单、配置灵活，处理简单、便于解释、正被广泛应用于全球的海洋地球物理调查中，如海域陆架区域地球物理勘查［3，4］，另外在新能源-天然气水合物的勘探［5，6］中也得到了应用。

单道地震接收电缆的好坏对于高质量的单道地震数据的采集有着重要的影响［7］，常见的接收电缆通常采用多个检波器组合的形式获取单道信号，因此检波器的灵敏度以及其组合方式等都会影响到采集数据的质量。另外，作业环境对于单道地震数据的采集也有着重要的影响［8］。

2 单道地震电缆设备以及数据比较分析

21 单道地震接收电缆性能比较

单道地震电缆内部的检波器通常是由压电材料制造，利用压电效应［9］将地震波引起的水压变化转变成电信号，然后通过电缆或者光缆等介质传输到船上的记录系统进行显示和存储。早期的海上地震勘探通常使用悬浮在固定或者缓慢漂流船只上的单个检波器来进行采集［10］。后来发展到特定排列方式组成的阵列拖缆拖在船尾，检波器被集成安置在充满油的柔韧塑料管中，保持与海水相适应的波阻抗，保证检波器与水的紧耦合。这一段包含检波器的塑料管就是本文所指的信号接收段，这种装置形式被设计成流线型来减小机器震动和水流引起的噪声。

本文所比较的两个电缆，一条是 Applied Acoustics Engineering公司生产的AAE20型单道接收电缆（以下简称为AAE电缆，图1左），其由20个检波器单元组成，间距015m，检波器段长度45m，其频率响应为145～7KHz（-3dB），总的灵敏度为-167 dB ref 1 v per。另外一条是Geo Marine Survey Systems公司生产的48检波器单元的Geo⁃Sense Mini⁃Streamers接收电缆（以下简称为GEO电缆，图1右），检波器间距04m，接收部分总长约20m，其检波器型号为AQ-2000，性能参数如下：灵敏度（@100Hz）为-201dB ref 1v per+/-15dB，频率响应为1Hz到10KHz（+/-2dB）。

图1 两种电缆实际形态（左边蓝色为AAE型电缆，右边红色为GEO型电缆）

Fig1 The samples oftwo cables（AAE Cable on the left，GEO Cable on the right）

22 数据采集及比较分析

本次野外数据采集作业水域水深在250m左右，统一采用GEO-6 KJ电火花震源系统作为激发源，选取合适且大小相同的激发能量，同时为了避免船尾螺旋桨引起的尾流以及环境噪声［11］对数据产生影响，采用如下装置参数：接收电缆释放长度选择为45m，炮检距（接收电缆与震源缆之间的距离）约为10m。另外，采集参数如下：采样频率（Sample frequency）4000Hz，延时（Delay）300ms，记录长度（Record Length）500ms，因此每道数据有2000个采样点

本文从野外实测数据中选取了一组进行比较，该组测线分别由上述两种不同的电缆进行数据采集，间距约300m，地质情况类似，且作业时测线方向一致，海况及作业环境相差不大，因此具有较强的可比性。

图2和图3分别显示GEO电缆和AAE电缆的原始实测数据，从图中可以看出两条电缆都有良好的穿透能力以及地层分辨能力，所采数据也都能满足常规单道地震数据的采集规范。

由于作业区域水相对较浅，AAE电缆检波器接收段较短且接收段内检波器之间距离较小，因此从图3的地震剖面可以看出该电缆接收到的数据在浅部有较高的分辨率。相比之下，GEO电缆检波器多且间距较大，其地震剖面信噪比高，尤其深部地层反射信号更加清晰。

图2 GEO电缆采集到的数据

Fig2 Data collected by GEO Cable

图3 AAE电缆采集到的数据

Fig3 Data collected by AAE Cable

图4和图5为从上述剖面中各抽取的一道数据，由于前置放大器设置不同造成两条电缆信号回波强度不一致，并不影响数据的信噪比及分辨率。

图4 GEO电缆抽取一道数据

Fig4 One channel extraction of GEO Cable

图5 AAE电缆抽取一道数据

Fig5 One channel extraction of AAE Cable

另外，分别对所抽取数据进行快速傅立叶变换（FFT），观测其频率域特征，图6和图7分别为上述两道数据的振幅谱，可以看到两条电缆能量主要集中在0-1000Hz之间，但是单从频率域并不能对两条电缆的分辨率、穿透能力以及信噪比等参数进行直观分析。

图6 GEO抽取数据的振幅谱

Fig6 The amplitude spectrogram of data extraction by GEO Cable

图7 AAE电缆抽取数据的振幅谱

Fig7 The amplitude spectrogram of data extraction by AAE Cable

因此，本文引入了S变换对所抽取数据进行了分析。传统的傅立叶变换虽然能对信号的整体性质进行分析，但却难获得信号的一些瞬时信息，特别是对于非平稳信号，而这些局部性质却是分析这些非平稳信号的关键。时频分析技术可以同时从时间域和频率域分析信号的局部特征，因此在傅立叶分析的基础上，发展了一系列的时频域信号分析论，如短时窗傅立叶变换、小波变换等。

S变换是由Stockwell（1996）等提出的，是连续小波变换的一种扩展，它的原理是基于一个高度和宽度随频率变化的高斯窗。函数的S变换定义如下：

南海地质研究（2014）

其中t表示时间，f表示频率，τ是一个控制高斯窗函数 在时间轴上位置的一个参数，e-2πft是一个相位因子，起到相位校正作用，这是较小波变换的一个优点［13］。另外，其离散形式可以利用快速傅立叶变换FFT的结果来计算，因此计算效率比较高。信号h（t）的离散形式h［kT］，k=0，1，…，N-1（T为采样间隔）的离散傅立叶变换表示如下：

南海地质研究（2014）

其中n=0，1…，N-1，则h（kT）的离散S变换可表示为：

南海地质研究（2014）

其中：j、m、n=0，1…，N-1。

图8和图9分别为GEO电缆和AAE电缆所抽取数据的S变换显示，从中可以看出两条电缆都能够采集到有效的地层反射信号，海底以及近海底地层反射信号最强，深部地层的反射信号也明显被记录到。另外，前述的AAE电缆浅层分辨率高，GEO电缆信噪比高也在图上有直观的反映。

3 结论及分析

本文采集了两条不同的单道地震接收电缆的数据，通过观察其剖面、频率域特征、时频域等特征，对上述两条电缆的特性、适用范围等进行了分析比较。

图8 GEO电缆抽取数据的S变换域显示

Fig8 The sample of S⁃transform domain of data extraction by GEO Cable

图9 AAE电缆抽取数据的S变换域显示

Fig9 The sample of S⁃transform domain of data extraction by AAE Cable

AAE电缆其检波器接收部分是由20个检波器组合形成，可以有效的压制环境噪声［14］。此外，由于其检波器间距只有015m，因此整个接收段部分长度仅为45m，相对较短，这种装置形式对于浅水区有明显的优势，在满足质量要求的前提下具有更高的分辨率。而本文所使用的GEO电缆其接收段较长，且检波器组合达到48个，除了达到压制噪声的目的外还可以有效采集到来自深水区以及深部地层的反射波，避免丢失回声数据，另外其具有较宽的频率响应以及良好的灵敏度，使其能采集到高信噪比的资料。

检波器组合具有低频响应［7］，因此会对高分辨率数据采集有一定影响，但是检波器组合又是高信噪比数据采集必不可少的，因此在野外作业时要根据作业环境、接收条件、工作目的等选择合适的接收电缆。

致谢

感谢参与数据采集的全体成员，中国地震局地球物理研究所吴庆举研究员提供了本文的S变换程序，温明明高工、吴衡高工、牟泽霖、林桂柱在本文的资料分析中给予悉心指导！

参考文献

［1］李军峰，肖都，孔广胜，等2004单道海上反射地震在海上物探工程中的应用［J］物探与化探，28（4）：365 368

［2］岳保静2010单道地震资料处理方法及应用［D］硕士学位论文，青岛：中国科学院海洋研究所

［3］D Van Rooij，B De Mol，V Huvenne，et al2003Seismic evidence of current⁃controlled sedimentation in the Belgica mound province，upper Porcupine slope，southwest of Ireland［J］Marine Geology，195（1-4）：31-53

［4］Nicolas Weber，Eric Chaumillon，Michel Tesson，et al2004Architecture and morphology of the outer segment of a mixed tide and wave⁃dominated⁃incised valley，revealed by HR seismic reflection profiling：the paleo⁃Charente River，France［J］Marine Geology，207（1-4）：17-38

［5］张光学，黄永祥，陈邦彦2003海域天然气水合物地震学［M］北京：海洋出版社

［6］吴志强，文丽，童思友，等2007海域天然气水合物的地震研究进展［J］地球物理学进展，22（1）：218-227

［7］付清锋2005地震检波器新技术发展方向［J］石油仪器，19（6）：1-4

［8］李丽青，梁蓓雯，徐华宁2007海上单道地震资料中多次波的衰减［J］石油物探，46（5）：457-464

［9］吴毅强，邓忠华1995压电器件［M］北京：电子工业出版社

［10］E J W Jones1999Marine Geophysics［M］John Wiley＆Sons

［11］刘建勋2007提高海上单道反射地震记录信噪比和分辨率的方法技术［J］物探化探计算技术，29（增刊）：116-122

［12］Stockwell R G，Manisinha，R P Lowe1996Localization of the complex spectrum：the S transform［J］IEEE TransSignal Process，44：998-1001

［13］Stefano Parolai2009Denoising of Seismograms Using the S Transform［J］Bulletin of the Seismological Society of America，99（1）：226-234

［14］陆基孟1993地震勘探原理（上、下册）［M］北京：石油大学出版社

Comparison Between Two Different Single-Channel Seismic Streamers

Feng Qiang Qiang1，2，Wen Ming Ming1，2，Wu Heng1，2

（1Guangzhou Marine Geological Survey，Guangzhou，510760；2Key laboratory of Marine Mineral Reasources，MLR，Guangzhou，510760）

Abstract：In this paper we collect the raw data from two different singlechannel reflection seismic streamers（AAE streamers and GEO streamers）and analyze those data comparatively from frequency domain and time-frequency domain by using fast Fourier transform（FFT）and StransformThe result shows that AAE streamers which have shorter receiving section and compact geophones can get higher shallow resolution and therefore are more suitable for shallow waterWhile GEO streamers are more suitable for deep water because of its longer receiving section，more geophones，and the ability to acquire high signal⁃to⁃noise ratio data with longer span between geophones

Key word：Single⁃channel reflection seismic；S⁃transform；Time⁃frequency domain；Geophone