数字信号处理的几个前沿题.docx
《数字信号处理的几个前沿题.docx》由会员分享，可在线阅读，更多相关《数字信号处理的几个前沿题.docx（23页珍藏版）》请在冰点文库上搜索。
数字信号处理的几个前沿题
第10章数字信号处理的几个前沿课题
前面介绍了数字信号处理的基本知识，本章我们将介绍时谱分析、小波变换、地震观测系统仿真与地面运动恢复等几个数字信号前沿课题，以便大家在实际工作中参考。
10.1时谱（倒谱）分析
时谱分析（Cepstrumanalysis）是一种非线性信号处理技术，它在语言、图像、和噪声处理领域中都有广泛的应用。
时谱可分为两类：
复时谱和功率时谱。
MATLAB信号处理工具箱提供复时谱分析的工具函数。
复时谱（Complexcepstrum）的定义为：
（10-1）
由上式可见，复时谱实际上是序列x（n）的Fourier变换取自然对数，再取Fourier逆变换，得到的复时谱仍然是一个序列。
也就是说，复时谱是x（n）从时间域至频率域、频率域至频率域、频率域至时间域的三次变换。
MATLAB信号处理工具箱函数cceps用于估计一个序列x的复时谱，调用格式为：
xhat=cceps（x）
式中，x为输入序列（实序列）；xhat为复时谱（复序列）。
MATLAB信号处理工具箱还提供了序列实时（倒）谱的计算程序rceps，调用格式为Y=rceps（x），其中x为实序列；y为实时谱，执行的操作为：
（10-2）
由此可知，我们不能从序列x的实时谱重构原始序列，因为实时谱是根据序列Fourier变换的幅值计算的，丢失了相位方面的信息。
但如果需要，可采用最小相位模式估计原始序列。
由于复时谱从复频谱计算得到，不损失相位信息，因此复时谱是可逆的，实时谱过程是不可逆的。
时谱分析技术广泛地应用于语言信号分析、同态滤波技术中。
这里举一个说明复时谱在具有回声信号测量中的应用。
【例10-1】设原信号是一个45Hz的正弦波，在传播过程中遇到障碍产生回声，回声振幅衰减为原信号的0.5，并与原信号有0.2s的延迟。
在某测点测到的信号是原信号和回声信号的叠加。
试使用复时谱分析该测点的信号。
%Samp10_1
t=0:
0.01:
1.49;%时间信号，采样间隔为0.01s
sig=sin（2*pi*45*t）;%原始信号
echo=0.5*[zeros（1,20）,sig（1:
length（t）-20）];
%在信号前面补加20个零，并使其振幅衰减一半作为回声信号
sigecho=sig+echo;
c=cceps（sigecho）;%求原始信号与回声信号叠加的复时谱
subplot（2,1,1）,plot（t,sigecho）;%绘制原始信号与回声的叠加信号
xlabel（'时间/s'）;axis（[01.5-11]）;gridon
title（'时间域信号'）
subplot（2,1,2）,plot（t,c）;%绘制其复时谱
xlabel（'时间/s'）;axis（[01.5-11]）;gridon
title（'时谱域信号'）
程序运行结果为图10-1。
可以看到，信号的复时谱在t=0.2秒处有一个峰值，这就是回声信号。
表示运用复时谱可以明确检测到回声信号。
在数字地震仪监测核爆破的工作中，如果采用两次核爆在几乎同一时间、同一地点发生，爆破产生的地震波经过地球内部传到同一地震台上的传播路径几乎一致，并且爆破源产生的波形相差不大，因此在地震台上表现为顺序到达的两个波形类似的波，若后面的波形被淹没，可以采用复时谱的方法进行检测，明确给出是否是两次爆破事件，如果是两次爆破事件，就可以明确给出两次事件的时间间隔。
图10-1正弦波与回声信号叠加的波形和时谱形状
10.2地震观测系统的仿真和地面运动的恢复
由前面的线性系统的讲解可知，一个系统可以用它的系统函数或脉冲响应来表示：
（10-3）
式中，x（t）为输入信号，对于地震观测系统来讲为地面运动；y（t）为系统的输出，对于地震观测系统来讲为地震记录；h（t）为系统的脉冲响应。
在频率域内，根据卷积定理，该式可以表示为：
（10-4）
式中，
为系统的传递函数，
和
分别为x（t）和y（t）的Fourier变换。
现在设想一个频带范围很宽的线性系统，比如宽频带地震仪，其系统函数为
;另一个频带较窄的系统，比如短周期地震仪，其系统函数为
，对于同样的输入
有：
，
（10-5）
式中，
为频带较窄的系统记录的频谱；
为频带较窄系统的传递函数。
比较前面两式可得：
（10-6）
将上式变换到时间域就得到频带较窄系统的输出y1（t）。
也就是说，如果知道了宽频带和窄频带系统的传递函数
和
，原则上可以从宽频带系统的输出推测出频带较窄系统的输出。
但如果我们知道窄频带系统的输出及其两种系统的传递函数，无法得到宽频带系统的输出。
这样就使得我们在记录某种信号时采用宽频带记录，然后仿真到其他各种窄频带的记录仪器上对信号进行分析。
如果已知地震仪的输出和地震仪的传递函数，我们可以求出地面运动为：
（10-7）
下面我们就举例给出将宽频带系统的输出仿真到窄频带输出及地面运动恢复的方法。
【例10-2】设计一个Butterworth模拟宽频带带通滤波器，设计指标为：
通带频率：
1~40Hz，低端阻带边界0.2Hz，高端阻带边界42.5Hz，通带波纹1dB，阻带衰减20dB。
再设计一个窄频带带通滤波器，设计指标为通带频率：
5~7Hz，低频段过渡带宽1.5Hz，高频段过渡带宽0.5Hz，通带波纹1dB，阻带衰减20dB。
假设一个信号
，其中f1=6Hz,f2=10Hz,f3=19Hz。
信号的采样频率为50Hz。
（1）模拟输入信号通过宽频带滤波器的输出信号并与原来信号进行比较。
（2）运用宽带滤波器的传递函数和输出得到输入信号，并与原输入信号进行比较（恢复地面运动）（3）模拟输入信号通过窄频带滤波器的输出信号并与原来信号进行比较。
（4）将在宽频带滤波器的模拟输出信号仿真到窄频带滤波器上并与（3）的结果进行比较。
%Samp10_2
Fs=100;dt=1/Fs;%给定模拟输入信号采样间隔
f1=6;f2=10;f3=19;%模拟输入信号的频率成分
N=500;%数据点数
t=[0:
N-1]*dt;%模拟输入信号的时间序列
x=sin（2*pi*f1*t）+0.5*cos（2*pi*f2*t）+0.5*sin（2*pi*f3*t）;%模拟输入信号
X=fft（x）;%得到输入信号的Fouier变换
%设计宽带Butterworth模拟滤波器
wp=[140]*2*pi;ws=[0.242.5]*2*pi;%通带和阻带截止频率
Rp=1;Rs=20;%通带波纹和阻带衰减
[Order,Wn]=buttord（wp,ws,Rp,Rs,'s'）;%求解Butterworth滤波器的最小阶数
w=linspace（0,Fs/2,N/2）*2*pi;%设置绘制幅频特性的频率
[b1,a1]=butter（Order,Wn,'s'）;%设计Butterworth带通滤波器
H=freqs（b1,a1,w）;%计算带通滤波器的频率响应
magH=abs（H）;phaH=unwrap（angle（H））;%求得振幅谱和相位谱
figure
（1）
subplot（3,1,1）,plot（w/（2*pi）,20*log10（magH））;%绘制幅频特性
xlabel（'频率/Hz'）
ylabel（'振幅/dB'）;
ylim（[-1000]）;xlim（[0,50]）;gridon
title（'宽带模拟带通滤波器'）;
%模拟输出
H=freqs（b1,a1,w）;%滤波器的幅频响应
H1=zeros（1,N）;
forii=1:
N/2%构建与Fourier变换得到前后数据共轭的形式
H1（ii）=H（ii）;
H1（N-ii）=conj（H1（ii））;
end
X=fft（x）;%将原信号进行Fourier变换
Y=zeros（1,N）;
forii=1:
N
Y（ii）=X（ii）.*H1（ii）;%按（10-4）式得到的输出
end
y=real（ifft（Y））;%将频率域的数据通过Fourier变换转换到时间域
subplot（3,1,2）,plot（t,x）,title（'输入信号'）%绘制输入信号
subplot（3,1,3）,plot（t,y）%绘制输出信号
title（'宽带滤波器的模拟输出信号'）
xlabel（'时间/s'）
%恢复地面运动
XX=zeros（1,N）;
forii=1:
N
if（abs（H1（ii））>1.0e-1）%为防止位于阻带内数据不稳定
XX（ii）=Y（ii）./H1（ii）;%恢复地面运动公式（10-7）
end
end
xx=real（ifft（XX））;%转换到时间域
figure
（2）
plot（t,x,t,xx,':
r'）
legend（'原始信号','恢复信号',1）%给出图例
title（'原始信号和恢复信号'）
xlabel（'时间/s'）
%窄带滤波器的设计
wp=[57]*2*pi;ws=[3.57.5]*2*pi;%窄带滤波器的通带和阻带边界频率
Rp=1;Rs=20;%窄带滤波器的通带波纹和阻带衰减
[Order,Wn]=buttord（wp,ws,Rp,Rs,'s'）;%计算窄带滤波器的最小阶数和截止频率
w=linspace（0,Fs/2,N/2）*2*pi;%设置绘制窄带滤波器幅频特性的频率
[b2,a2]=butter（Order,Wn,'s'）;%设计Butterworth窄带滤波器
H=freqs（b2,a2,w）;%计算Butterworth窄带滤波器的频率响应
magH=abs（H）;phaH=unwrap（angle（H））;%计算振幅谱和相位谱
figure（3）
subplot（3,1,1）,plot（w/（2*pi）,20*log10（magH））;%绘制幅频特性
xlabel（'频率/Hz'）,ylabel（'振幅/dB'）;ylim（[-1000]）;xlim（[0,50]）
title（'窄带模拟滤波器'）;
gridon
%模拟输出
H=freqs（b2,a2,w）;
H2=zeros（1,N）;
forii=1:
N/2
H2（ii）=H（ii）;
H2（N-ii）=conj（H2（ii））;
end
X=fft（x）;
Y1=zeros（1,N）;
forii=1:
N
Y1（ii）=X（ii）.*H2（ii）;
end
y1=real（ifft（Y1））;
subplot（3,1,2）,plot（t,x）,title（'输入信号'）%绘制输入信号
subplot（3,1,3）,plot（t,y1）%绘制输出信号
title（'窄带滤波器的模拟输出'）
xlabel（'时间/s'）
%运用宽频带仪器的输出仿真得到窄带仪器上的输出
figure（4）
XY=zeros（1,N）;
forii=1:
N
if（abs（H1（ii））>1.0e-2）
XY（ii）=Y（ii）*H2（ii）/H1（ii）;%公式（10-6）
end
end
xy=real（ifft（XY））;%转换到时间域
plot（t,y1,t,xy,':
r'）
legend（'窄带输出','仿真输出',1）
xlabel（'时间/s'）
程序的输出为图10-2~5。
可见运用宽频带的仪器响应和其输出信号，可以较为精确地恢复地面运动（图10-3）；将宽带滤波器的输入仿真到窄带滤波器上，得到和窄带滤波器一样的波形（图10-5），验证了上述方法的正确性。
图10-2例10-2中宽频带仪器的幅频响应、输入信号和输出信号
图10-3例10-2中利用宽频带仪器的输出和仪器的传递函数得到输入信号并与原输入信号进行比较
图10-4例10-2中窄带仪器的幅频响应、输入信号和输出信号
图10-5例10-2中运用宽带仪器的输出、宽带和窄带仪器的传递函数得到的窄带仪器的输出（仿真），并与原窄带仪器的输出进行比较
下面我们用真实的地震数据和数字滤波器方法来揭示地面运动恢复和向窄带地震仪上的仿真。
【例10-3】以首都圈地震记录文件1653.VBB的hns台站的上下向记录为例，运用椭圆滤波器来揭示地面运动恢复和仿真的概念。
宽带仪器的通带边界频率为[0.00124.8]Hz，阻带边界频率为[0.0000124.9]Hz，通带波纹为1dB，阻带衰减为50dB。
相当于短周期地震仪的窄带仪器的阻带边界频率为[0.014.5]Hz,通带边界频率为[0.13.8]Hz，通带波纹为1dB，阻带衰减为20dB；相当于长周期地震仪的窄带仪器用低通滤波器来表示，阻带边界频率为0.1Hz,通带边界频率为0.02Hz，通带波纹为1dB，阻带衰减为30dB。
（1）采用宽频带仪器传递函数和输出信号，恢复输入信号（地面运动）；
（2）将宽频带仪器的输出仿真到短周期窄带仪器上；并与窄带仪器的输出进行比较；（3）将宽频带仪器的输出仿真到长周期窄带仪器上；并与窄带仪器的输出进行比较；
程序如下：
%Samp10_3
loadhns.dat;%假设的地震仪记录之前的地面运动
Xt=hns;%把数据赋值给变量
Fs=50;%设定采样率单位（HZ）
dt=1/Fs;%求采样间隔单位（s）
N=length（Xt）;%得到序列的长度
t=[0:
N-1]*dt;%时间序列
%设计一个椭圆滤波器，假定为宽频带地震仪
ws=[0.0000124.9]*2/Fs;wp=[0.00124.8]*2/Fs;%通带和阻带边界频率（归一化频率）
Rp=1;Rs=50;Nn=512;%通带波纹和阻带衰减以及绘制频率特性的数据点数
[Order,Wn]=ellipord（wp,ws,Rp,Rs）;%求取数字滤波器的最小阶数和归一化截止频率
[b,a]=ellip（Order,Rp,Rs,Wn）;%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器
figure
（1）
[H,f]=freqz（b,a,Nn,Fs）;%按传递函数系数、数据点数和采样频率求得滤波器的频率特性
subplot（2,1,1）,plot（f,20*log10（abs（H）））
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,1,2）,plot（f,180/pi*unwrap（angle（H）））
xlabel（'频率/Hz'）;ylabel（'相位/^o'）;gridon;
y=filtfilt（b,a,Xt）;%在宽带滤波器上的输出
figure
（2）
subplot（2,1,1）,plot（t,Xt）
xlabel（'时间/s'）,title（'输入信号'）;
ylabel（'振幅'）;
subplot（2,1,2）,plot（t,y）
xlabel（'时间/s'）,title（'输出信号'）;
ylabel（'振幅'）;
%已知宽频带地震仪的频率特性,恢复地面运动
[H,f]=freqz（b,a,N,Fs,'whole'）;%得到地震仪的特性
Y=fft（y）;
XX=zeros（1,N）;
forii=1:
N%按（10-7）式得到地面运动的频率域表示
if（H（ii）>1.0e-4）
XX（ii）=Y（ii）./H（ii）;
end
end
disp=real（ifft（XX））;%得到地面运动
figure（3）,plot（t,Xt,t,disp,'r'）
legend（'原始信号','恢复地面运动',1）
xlabel（'时间/s'）
ylabel（'振幅'）;
%仿真到短周期地震仪上，短周期地震仪用一个窄带椭圆滤波器来表示
ws=[0.014.5]*2/Fs;wp=[0.13.8]*2/Fs;%通带和阻带边界频率（归一化频率）
Rp=1;Rs=20;Nn=512;%通带波纹和阻带衰减以及绘制频率特性的数据点数
[order,Wn]=ellipord（wp,ws,Rp,Rs）;%求取数字滤波器的最小阶数和归一化截止频率
[b,a]=ellip（order,Rp,Rs,Wn）;%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器
figure（4）
[H1,f]=freqz（b,a,Nn,Fs）;%按传递函数系数、数据点数和采样频率求得滤波器的频率特性
subplot（2,1,1）,plot（f,20*log10（abs（H1）））
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,1,2）,plot（f,180/pi*unwrap（angle（H1）））
xlabel（'频率/Hz'）;ylabel（'相位/^o'）;gridon;
figure（5）
y1=filtfilt（b,a,Xt）;%在窄带滤波器上的输出
[H1,f]=freqz（b,a,N,Fs,'whole'）;%得到地震仪的特性
XX1=zeros（1,N）;
forii=1:
N%按（10-6）式得到仿真结果
if（abs（H（ii））>1.0e-4）
XX1（ii）=Y（ii）.*H1（ii）/H（ii）;
end
end
x1=ifft（XX1）;
plot（t,y1,t,real（x1）,'r'）%绘制输入信号
legend（'实际输出','仿真输出',1）
xlabel（'时间/s'）;
ylabel（'振幅'）;
gridon;
%仿真到长周期地震仪上，长周期地震仪用一个窄带椭圆滤波器来表示
ws=0.1*2/Fs;wp=0.02*2/Fs;%通带和阻带边界频率（归一化频率）
Rp=1;Rs=30;Nn=512;%通带波纹和阻带衰减以及绘制频率特性的数据点数
[Order,Wn]=ellipord（wp,ws,Rp,Rs）;%求取数字滤波器的最小阶数和归一化截止频率
[b,a]=ellip（Order,Rp,Rs,Wn）;%按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器
figure（6）
y1=filtfilt（b,a,Xt）;%在窄带滤波器上的输出
[H1,f]=freqz（b,a,N,Fs,'whole'）;%得到地震仪的特性
XX1=zeros（1,N）;
forii=1:
N
if（abs（H（ii））>1.0e-4）%为了防止H值太小将该频率的信号放大
XX1（ii）=Y（ii）.*H1（ii）./H（ii）;%按（10-6）式得到仿真结果
end
end
x1=ifft（XX1）;
plot（t,y1,t,real（x1）,'r:
'）%绘制输入信号
legend（'实际输出','仿真输出',1）
xlabel（'时间/s'）;
ylabel（'振幅'）;
gridon;
图10-6例10-3中宽带仪器的频率特性：
上图：
幅频响应；下图：
相频响应
图10-7例10-3中宽带仪器输出（下图）和原输入信号（上图）的比较
图10-8例10-3中运用宽带仪器传递函数和宽带仪器的输出信号恢复的地面运动和原输入信号的比较
图10-9例10-3中短周期窄带仪器的频率特性：
上图：
幅频响应；下图：
相频响应
图10-10例10-3中运用宽带和短周期窄带仪器传递函数及宽带仪器输出仿真到窄带仪器上并与窄带仪器的实际输出进行比较
图10-11例10-3中运用宽带和长周期窄带仪器传递函数及宽带仪器输出仿真到窄带仪器上并与窄带仪器的实际输出进行比较
图10-6为宽带仪器的频率响应，我们设计了较宽的频带，几乎包含了输入信号的所有频率。
输入信号经过宽带仪器输出后波形基本未变，但幅值发生了一些改变（图10-7）。
运用宽带仪器的输出及传递函数恢复的地面运动与原输入信号相差很小（图10-8）。
图10-9为相对的短周期窄带仪器的频率特性，运用宽带仪器输出信号、宽带和窄带仪器的传递函数仿真到窄带仪器上与原信号经窄带仪器的滤波结果相差很小（图10-10），由于是短周期仪器，因此仿真结果中短周期体波成分较为丰富，几乎看不到面波成分。
运用宽带仪器输出信号、宽带和窄带长周期仪器的传递函数仿真到窄带仪器上与原信号经窄带仪器的滤波结果也相差不大（图10-11），由于是长周期仪器，因此仿真结果中短周期体波成分很小，长周期面波成分较为丰富。
10.3小波分析举例
小波变换最早是1984年研究地震资料而使用的。
由于小波分析具有能够根据分析对象自动调整有关参数的“自适应性”和能够根据观测对象自动“调焦”的特性而广泛应用于地震震相识别、地震波形反演等研究领域。
一个地震时间序列f（t）是由不同强度的大小地震构成的，每个时刻t的强度可以用下式表示：
（10-8）
其中
为奇异性指数或标度指数。
奇异性是指
时
的数值（特别是
）决定了
的速度。
标度性是指若变量t变成
，则
和
是自相似的，即：
（10-9）
（10-9）式常称为标度律。
自相似的标度律是多种尺度（即无特征尺度）现象的特征。
由于中国地震资料的时间序列含有多种不同尺度，也即包含有很多层次。
虽然
是杂乱无章的，但小波变换好比显微镜，调节其中的放大倍数就可以清楚地看出
在各个层次上的变化趋势。
由于未对地震资料作谱分析，所以这是地震资料的较为真实的物理分析。
10.3.1小波变换与突变
一个时间序列信号
的Fourier变换，是从频率轴（或波数域）来分析该信号是由哪些频率（或波数）的波组成的。
地震资料
组成的序列并不是平稳的，而是随时间不断变化的，此时Fourier变换就不能提取
中的奇异性和突变点的信息，它只是将这些信息铺开到整个频率轴上。
突变点对于很多问题是很重要的，例如，全球气候变暖从何时开始，地震活动性何时开始表现活跃等。
但是小波变换Tg是将
分解成具有局部特性的小波
:
（10-10）
其中小波
是将具有局部特性的小波函数
通过平移
和放大（放大倍数为1/a）而构成的。
参数a具有时间的量纲，也称为小波尺度。
设
为具有低通性质的平滑函数，以它的一阶、二阶导数作为小波对f（t）作小波变换，可以证明:
（10-11）
（10-12）
式中，
表示卷积。
也就是说，用
、
对f（t）作小波变换分别相当于f（t）被
平滑后在对t求一阶或二阶导数。
因此，对某一固定a值，f（t）
的拐点即是
的极值点，又是
的过零点，由此可检测信号的急剧变化之处。
本研究中所采用低通性质的平滑函数为高斯函数，根据它的一阶导数、二阶导数（墨西哥草帽小波函数）作为小波基函数进行突变点分析。
各小波函数的表达式为：
（10-13）
之所以选择此函数是因为它对称、可微、可积，时频两域都是高斯型且呈平方型指数衰减特性，在时频两域均具有很好的局域性。
该函数的MATLAB程序如下：
function[x,y,z]=gdttest（a,t0,si）
%wavelettransformMexscroHat
t=（1:
1:
si）;%给出时间序列
x=（1-（（t-t0）/a）.^2）.*exp（-（（t-t0）/a）.^2/2）/a;%二阶导数
y=exp（-（（t-t0）/a）.^2/2）/a;
升级会员 