震波观测系统MATLAB仿真
课程名称 数字信号处理
实验项目 题目6 震波观测系统MATLAB仿真
指导教师
学 院 光电信息通信工程 _
专 业 电子信息工程
班级学号
学生姓名
课设时间 201112 28201215
数字信号处理课程设计务书
题目6
震波观测系统MATLAB仿真
容
掌握震波观测系统数字信号处理方法实现宽频带系统输出仿真窄频带输出面运动恢复
设计
求
求
某震台站记录震观测文件例选择合适滤波器揭示面运动恢复仿真概念
步骤
1读取震波观测文件数做出时域频域图形设计包含频率成分宽频带滤波器假定宽频带震仪恢复面运动绘出滤波器频率特性面运动时域图
2已知短周期窄带仪器阻带边界频率[001 45]Hz通带边界频率[01 38]Hz通带波纹1dB阻带衰减20dB 宽频带仪器输出仿真短周期窄带仪器窄带仪器输出进行较（画图）绘出窄带仪器频谱图
3长周期震仪窄带仪器低通滤波器表示阻带边界频率01Hz通带边界频率002Hz通带波纹1dB阻带衰减30dB宽频带仪器输出仿真长周期窄带仪器窄带仪器输出较步骤2作图
仪
器设备
1计算机1台安装MATLAB软件
参
考文献
[美]数字信号处理——MATLAB[M]西安：西安交通学出版社2002
课程设计进度计划（起止时间工作容）
课程设计安排6题目中题目整课程设计24学时分15周安排具体进度：
4学时 复题目相关知识掌握实现原理
12学时 MATLAB语言实现题目求
4学时 进步完善功现场检查答辩
4学时 完成课程设计报告
课程设计开始日期
20111226
课程设计完成日期
201216
课程设计实验室名称
信号信息处理实验室
点
实验楼3603605
资料载址
目录
摘 4
正文 4
目 4
二原理 4
三求 5
四步骤 5
五程序实现 6
实验结果 12
六体会 15
参考文献 15
摘
文目实现震波观测系统MATLAB仿真线性系统y(t)h(t)*x(t)x(t)面运动h(t)系统击响应y(t)系统输出根卷积定理Y（ω）H(ω)X(ω)震波观测文件数y（t）设计宽频带滤波器h（t）恢复面运动x（t）短周期震仪系统函数H1（w）输入面运动x（t）Y1(ω)H1(ω)X（ω）推导出Y1（w）H1(w)Y(w)H(w)Y1（w）作ifft实现宽频带仪器短周期窄带仪器仿真样长周期震仪系统函数H2（w）Y2（w）H1(w)Y(w)H(w)然Y2（w）作ifft实现仿真椭圆滤波器巴特沃斯滤波器切雪夫滤波器设计简单滤波器指标没问题调相应函数实现仿真结果请参考文正文部分
正文
目
运学数字信号处理基知识掌握震波观测系统数字信号处理方法实现宽频带系统输出仿真窄频带输出面运动恢复
二原理
线性系统系统函数脉响应表示
y(t)h(t)*x(t) ①
式中x(t)输入信号相震观测系统面运动y(t)系统输出相震观测系统震记录h(t)系统击响应频率域根卷积定理该式表示
Y（ω）H(ω)X(ω) ②
式中H(ω)系统传递函数 X(ω) Y（ω）x(t)y(t)傅里叶变换
设想频带范围宽线性系统宽频带震仪系统函数H(ω)频带较窄系统短周期震仪系统函数H1(ω)样输入X(ω)
Y(ω)H(ω)X(ω) Y1(ω)H1(ω)X（ω） ③
式中Y1(ω)频带较窄系统记录频谱H1(ω频带较窄系统传递函数式③
H1(ω) Y(ω)
Y1(ω) H(ω) ④
式变换时间域频带较窄系统输出y1(t)说果知道宽频带窄频带系统传递函数H(ω) H1(ω)原宽频带系统输出推测出窄频带系统输出果知道窄频带系统输出两种系统传递函数法宽频带系统输出样记录某种信号时采宽频带记录然仿真种窄频带记录仪器信号进行分析
果已知震仪输出震仪传递函数求出面运动
X（ω） Y(ω) H(ω) ⑤
三求
某震台站记录震观测文件例选择合适滤波器揭示面运动恢复仿真概念
四步骤
1 读取震波观测文件数做出时域频域图形设计包含频率成分宽频带滤波器假定宽频带震仪恢复面运动绘出滤波器频率特性面运动时域图
2 已知短周期窄带仪器阻带边界频率[001 45]Hz通带边界频率[01 38]Hz通带波纹1dB阻带衰减20dB宽频带仪器输出仿真短周期窄带仪器窄带仪器输出进行较（画图）绘出窄带仪器频谱图长周期震仪窄带仪器低通滤波器表示阻带边界频率01Hz通带边界频率002Hz通带波纹1dB阻带衰减30dB宽频带仪器输出仿真长周期窄带仪器窄带仪器输出较步骤2作图
五程序实现
close allclear allclc
load hnsdat 读取数序列
Xthns 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(1)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
xlabel(＇时间s＇)title(＇时间域＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf))绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
设计切雪夫1型宽频带滤波器假定宽频带震仪
ws[000001 250]*2Fs 阻带边界频率（化频率）
wp[0001 250]*2Fs 通带边界频率（化频率）
Rp1Rs20Nn513 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]cheb1ord(wpwsRpRs)求取数字滤波器阶数化截止频率
[ba]cheby1(OrderRpWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(2)
[Hf]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
y1filtfilt(baXt)
subplot(211)plot(f20*log10(abs(H))) 画出宽带滤波器幅频特性
xlabel(＇＼lambda＇)ylabel(＇A(＼lambda)db＇)
title(＇宽频带滤波器幅频特性＇)grid on
subplot(212)plot(fangle(H)) 画出宽带滤波器相频特性
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)title(＇宽频带滤波器相频特性＇)grid on
已知宽频带震仪频率特性恢复面运动
[Hf]freqz(baNFs＇whole＇) 震仪特性
Xfzeros(1N)
for i1N
if (H(i)>10e4)
Xf(i)Yf(i)H(i) 面运动频率域表示
end
end
figure(3)
xtreal(ifft(Xf)) 面运动
subplot(211)
plot(txt＇r＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
title(＇面运动时域图＇)
grid on
subplot(212)
plot(tXt＇g＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
title(＇原始信号＇)
grid on
设计椭圆宽带滤波器假定宽频带震仪
ws[000001 250]*2Fswp[0001 250]*2Fs 通带阻带边界频率（化频率）
Rp1Rs50Nn512 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]ellipord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]ellip(OrderRpRsWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(4)
[Hf]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
subplot(211)plot(f20*log10(abs(H)))
xlabel(＇频率Hz＇)ylabel(＇振幅dB＇)grid on
subplot(212)plot(f180pi*unwrap(angle(H)))
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)grid on
yfiltfilt(baXt) 宽带滤波器输出
figure(5)
subplot(211)plot(tXt)
xlabel(＇时间s＇)title(＇输入信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
figure(6)
subplot(211)plot(ty1＇g＇)
xlabel(＇时间s＇)title(＇切雪夫1型宽频带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty＇r＇)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
仿真长周期震仪长周期震仪巴特沃思滤波器表示
ws01*2Fswp002*2Fs 通带阻带边界频率（化频率）
Rp1Rs30Nn512 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]buttord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]butter(OrderWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(7)
[H2f]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
subplot(211)plot(f20*log10(abs(H2)))
xlabel(＇＼lambda＇)ylabel(＇A(＼lambda)db＇)title(＇长周期窄带滤波器幅频特性＇)grid on
subplot(212)plot(fangle(H2))
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)title(＇长周期窄带滤波器相频特性＇)grid on
figure(8)
y2filtfilt(baXt) 窄带滤波器输出
[H2f]freqz(baNFs＇whole＇) 震仪特性
Yf2zeros(1N)
for i1N
if (abs(H2(i))>10e4) 防止H值太该频率信号放
Yf2(i)Yf(i)*H2(i)H(i) 仿真结果
end
end
x2ifft(Yf2)
subplot(211)
plot(ty2＇g＇) 绘制实际输出信号
xlabel(＇时间s＇)
ylabel(＇振幅＇)
title(＇长周期震仪实际输出＇)
grid on
subplot(212)
plot(treal(x2)＇r＇) 绘制仿真输出信号
title(＇长周期震仪仿真输出＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
仿真长周期震仪长周期震仪窄带椭圆滤波器表示
ws01*2Fswp002*2Fs 通带阻带边界频率（化频率）
Rp1Rs30Nn512 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]ellipord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]ellip(OrderRpRsWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(9)
y1filtfilt(baXt) 窄带滤波器输出
[H1f]freqz(baNFs＇whole＇) 震仪特性
XX1zeros(1N)
for ii1N
if (abs(H1(ii))>10e4) 防止H值太该频率信号放
XX1(ii)Yf(ii)*H1(ii)H(ii) 仿真结果
end
end
x1ifft(XX1)
subplot(121)
plot(ty1)
title(＇实际输出＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
subplot(122)
plot(treal(x1))
title(＇仿真输出＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
figure(10)
subplot(211)plot(ty2＇g＇)
xlabel(＇时间s＇)title(＇巴特沃思滤波器滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty1＇r＇)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
仿真短周期震仪短周期震仪窄带椭圆滤波器表示
ws[001 45]*2Fswp[01 38]*2Fs 通带阻带边界频率（化频率）
Rp1Rs20Nn512 通带波纹阻带衰减绘制频率特性数点数
[orderWn]ellipord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]ellip(orderRpRsWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(11)
[H1f]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
subplot(211)plot(f20*log10(abs(H1)))
xlabel(＇频率Hz＇)ylabel(＇振幅dB＇)grid on
subplot(212)plot(f180pi*unwrap(angle(H1)))
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)grid on
figure(12)
y1filtfilt(baXt) 窄带滤波器输出
[H1f]freqz(baNFs＇whole＇) 震仪特性
XX1zeros(1N)
for ii1N 仿真结果
if (abs(H1(ii))>10e4)
XX1(ii)Yf(ii)*H1(ii)H(ii)
end
end
x1ifft(XX1)
plot(ty1treal(x1)＇r＇) 绘制输入信号
legend(＇实际输出＇＇仿真输出＇1)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
仿真短周期震仪短周期震仪切雪夫2型滤波器表示
ws[001 45]*2Fswp[01 38]*2Fs
Rp1Rs20Nn512
[OrderWn]cheb2ord(wpwsRpRs)求取数字滤波器阶数化截止频率
[ba]cheby2(OrderRpWn)阶数截止频率通带波纹阻带衰减设计滤波器
figure(13)
[Hf]freqz(baNnFs)传递函数系数数点数采样频率求滤波器频率特性
y3filtfilt(baXt)
subplot(211)plot(f20*log10(abs(H)))画出宽带滤波器幅频特性
xlabel(＇＼lambda＇)ylabel(＇A(＼lambda)db＇)
title(＇宽频带滤波器幅频特性＇)
grid on
subplot(212)plot(fangle(H)) 画出宽带滤波器相频特性
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)
title(＇宽频带滤波器相频特性＇)
grid on
figure(14)
subplot(211)plot(ty3＇g＇)
xlabel(＇时间s＇)title(＇切雪夫2型滤波器滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty1＇r＇)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
close allclear allclc
load hns1dat 读取数序列
Xthns1 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(1)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
title(＇P波＇)
xlabel(＇时间s＇)title(＇时间域＇)
title(＇P波＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf)) 绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
load hns2dat 读取数序列
Xthns2 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(2)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
title(＇S波＇)
xlabel(＇时间s＇)title(＇时间域＇)
title(＇S波＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf)) 绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
load hns3dat 读取数序列
Xthns3 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(3)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
title(＇面波＇)
xlabel(＇时间s＇)title(＇时间域＇)
title(＇面波＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf))绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
实验结果
图1 切雪夫1型宽频带滤波器椭圆宽带滤波器输出信号
面运动时域原始信号
图2 输入信号输出信号
图3 宽频带振幅相位
图4 短周期窄带振幅相位
图5 实际输出仿真输出
图6 巴特沃夫长周期实际输出仿真输出
图7 震波面波P波S波幅频图
图8 长周期窄带滤波器幅频特性
长周期窄带滤波器相频特性
六实验体会
通次实验进步复数字信号处理关滤波器基础解理实际身边处处数字信号处理相关知识应语音识采集处理震波观测直观证实数字信号处理门课程重性
次实验中实际操作中加强实践力巩固数字信号处理理知识培养解决实际问题力设计程中提高思考力动手力学理知识时明白应实际
次课程设计认识足认识学基础知识究竟运什领域运老师学耐心指导发现选择巴特沃斯切雪夫滤波器问题修改调试终应效果理实践相结合优势处受益匪浅
参考文献
[1]焦瑞莉 罗倩 汪毓铎 顾奕数字信号处理[M]机械工艺出版社pp184195
[2]httplgbougzcomcnhtmlautooldprotel99yyong1htm
袁宇波MatlabProtel设计微机保护中Butterworth模拟低通滤波器
[3]庆禹电力系统数字光电量测系统原理技术[J]电网技术200125(4) pp15
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容，不代表《香当网》观点或立场，请自行判断内容的真实性和可靠性！
该内容是文档的文本内容，更好的格式请下载文档
震波观测系统MATLAB仿真
课程名称 数字信号处理
实验项目 题目6 震波观测系统MATLAB仿真
指导教师
学 院 光电信息通信工程 _
专 业 电子信息工程
班级学号
学生姓名
课设时间 201112 28201215
数字信号处理课程设计务书
题目6
震波观测系统MATLAB仿真
容
掌握震波观测系统数字信号处理方法实现宽频带系统输出仿真窄频带输出面运动恢复
设计
求
求
某震台站记录震观测文件例选择合适滤波器揭示面运动恢复仿真概念
步骤
1读取震波观测文件数做出时域频域图形设计包含频率成分宽频带滤波器假定宽频带震仪恢复面运动绘出滤波器频率特性面运动时域图
2已知短周期窄带仪器阻带边界频率[001 45]Hz通带边界频率[01 38]Hz通带波纹1dB阻带衰减20dB 宽频带仪器输出仿真短周期窄带仪器窄带仪器输出进行较（画图）绘出窄带仪器频谱图
3长周期震仪窄带仪器低通滤波器表示阻带边界频率01Hz通带边界频率002Hz通带波纹1dB阻带衰减30dB宽频带仪器输出仿真长周期窄带仪器窄带仪器输出较步骤2作图
仪
器设备
1计算机1台安装MATLAB软件
参
考文献
[美]数字信号处理——MATLAB[M]西安：西安交通学出版社2002
课程设计进度计划（起止时间工作容）
课程设计安排6题目中题目整课程设计24学时分15周安排具体进度：
4学时 复题目相关知识掌握实现原理
12学时 MATLAB语言实现题目求
4学时 进步完善功现场检查答辩
4学时 完成课程设计报告
课程设计开始日期
20111226
课程设计完成日期
201216
课程设计实验室名称
信号信息处理实验室
点
实验楼3603605
资料载址
目录
摘 4
正文 4
目 4
二原理 4
三求 5
四步骤 5
五程序实现 6
实验结果 12
六体会 15
参考文献 15
摘
文目实现震波观测系统MATLAB仿真线性系统y(t)h(t)*x(t)x(t)面运动h(t)系统击响应y(t)系统输出根卷积定理Y（ω）H(ω)X(ω)震波观测文件数y（t）设计宽频带滤波器h（t）恢复面运动x（t）短周期震仪系统函数H1（w）输入面运动x（t）Y1(ω)H1(ω)X（ω）推导出Y1（w）H1(w)Y(w)H(w)Y1（w）作ifft实现宽频带仪器短周期窄带仪器仿真样长周期震仪系统函数H2（w）Y2（w）H1(w)Y(w)H(w)然Y2（w）作ifft实现仿真椭圆滤波器巴特沃斯滤波器切雪夫滤波器设计简单滤波器指标没问题调相应函数实现仿真结果请参考文正文部分
正文
目
运学数字信号处理基知识掌握震波观测系统数字信号处理方法实现宽频带系统输出仿真窄频带输出面运动恢复
二原理
线性系统系统函数脉响应表示
y(t)h(t)*x(t) ①
式中x(t)输入信号相震观测系统面运动y(t)系统输出相震观测系统震记录h(t)系统击响应频率域根卷积定理该式表示
Y（ω）H(ω)X(ω) ②
式中H(ω)系统传递函数 X(ω) Y（ω）x(t)y(t)傅里叶变换
设想频带范围宽线性系统宽频带震仪系统函数H(ω)频带较窄系统短周期震仪系统函数H1(ω)样输入X(ω)
Y(ω)H(ω)X(ω) Y1(ω)H1(ω)X（ω） ③
式中Y1(ω)频带较窄系统记录频谱H1(ω频带较窄系统传递函数式③
H1(ω) Y(ω)
Y1(ω) H(ω) ④
式变换时间域频带较窄系统输出y1(t)说果知道宽频带窄频带系统传递函数H(ω) H1(ω)原宽频带系统输出推测出窄频带系统输出果知道窄频带系统输出两种系统传递函数法宽频带系统输出样记录某种信号时采宽频带记录然仿真种窄频带记录仪器信号进行分析
果已知震仪输出震仪传递函数求出面运动
X（ω） Y(ω) H(ω) ⑤
三求
某震台站记录震观测文件例选择合适滤波器揭示面运动恢复仿真概念
四步骤
1 读取震波观测文件数做出时域频域图形设计包含频率成分宽频带滤波器假定宽频带震仪恢复面运动绘出滤波器频率特性面运动时域图
2 已知短周期窄带仪器阻带边界频率[001 45]Hz通带边界频率[01 38]Hz通带波纹1dB阻带衰减20dB宽频带仪器输出仿真短周期窄带仪器窄带仪器输出进行较（画图）绘出窄带仪器频谱图长周期震仪窄带仪器低通滤波器表示阻带边界频率01Hz通带边界频率002Hz通带波纹1dB阻带衰减30dB宽频带仪器输出仿真长周期窄带仪器窄带仪器输出较步骤2作图
五程序实现
close allclear allclc
load hnsdat 读取数序列
Xthns 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(1)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
xlabel(＇时间s＇)title(＇时间域＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf))绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
设计切雪夫1型宽频带滤波器假定宽频带震仪
ws[000001 250]*2Fs 阻带边界频率（化频率）
wp[0001 250]*2Fs 通带边界频率（化频率）
Rp1Rs20Nn513 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]cheb1ord(wpwsRpRs)求取数字滤波器阶数化截止频率
[ba]cheby1(OrderRpWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(2)
[Hf]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
y1filtfilt(baXt)
subplot(211)plot(f20*log10(abs(H))) 画出宽带滤波器幅频特性
xlabel(＇＼lambda＇)ylabel(＇A(＼lambda)db＇)
title(＇宽频带滤波器幅频特性＇)grid on
subplot(212)plot(fangle(H)) 画出宽带滤波器相频特性
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)title(＇宽频带滤波器相频特性＇)grid on
已知宽频带震仪频率特性恢复面运动
[Hf]freqz(baNFs＇whole＇) 震仪特性
Xfzeros(1N)
for i1N
if (H(i)>10e4)
Xf(i)Yf(i)H(i) 面运动频率域表示
end
end
figure(3)
xtreal(ifft(Xf)) 面运动
subplot(211)
plot(txt＇r＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
title(＇面运动时域图＇)
grid on
subplot(212)
plot(tXt＇g＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
title(＇原始信号＇)
grid on
设计椭圆宽带滤波器假定宽频带震仪
ws[000001 250]*2Fswp[0001 250]*2Fs 通带阻带边界频率（化频率）
Rp1Rs50Nn512 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]ellipord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]ellip(OrderRpRsWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(4)
[Hf]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
subplot(211)plot(f20*log10(abs(H)))
xlabel(＇频率Hz＇)ylabel(＇振幅dB＇)grid on
subplot(212)plot(f180pi*unwrap(angle(H)))
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)grid on
yfiltfilt(baXt) 宽带滤波器输出
figure(5)
subplot(211)plot(tXt)
xlabel(＇时间s＇)title(＇输入信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
figure(6)
subplot(211)plot(ty1＇g＇)
xlabel(＇时间s＇)title(＇切雪夫1型宽频带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty＇r＇)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
仿真长周期震仪长周期震仪巴特沃思滤波器表示
ws01*2Fswp002*2Fs 通带阻带边界频率（化频率）
Rp1Rs30Nn512 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]buttord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]butter(OrderWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(7)
[H2f]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
subplot(211)plot(f20*log10(abs(H2)))
xlabel(＇＼lambda＇)ylabel(＇A(＼lambda)db＇)title(＇长周期窄带滤波器幅频特性＇)grid on
subplot(212)plot(fangle(H2))
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)title(＇长周期窄带滤波器相频特性＇)grid on
figure(8)
y2filtfilt(baXt) 窄带滤波器输出
[H2f]freqz(baNFs＇whole＇) 震仪特性
Yf2zeros(1N)
for i1N
if (abs(H2(i))>10e4) 防止H值太该频率信号放
Yf2(i)Yf(i)*H2(i)H(i) 仿真结果
end
end
x2ifft(Yf2)
subplot(211)
plot(ty2＇g＇) 绘制实际输出信号
xlabel(＇时间s＇)
ylabel(＇振幅＇)
title(＇长周期震仪实际输出＇)
grid on
subplot(212)
plot(treal(x2)＇r＇) 绘制仿真输出信号
title(＇长周期震仪仿真输出＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
仿真长周期震仪长周期震仪窄带椭圆滤波器表示
ws01*2Fswp002*2Fs 通带阻带边界频率（化频率）
Rp1Rs30Nn512 通带波纹阻带衰减绘制频率特性数点数
[OrderWn]ellipord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]ellip(OrderRpRsWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(9)
y1filtfilt(baXt) 窄带滤波器输出
[H1f]freqz(baNFs＇whole＇) 震仪特性
XX1zeros(1N)
for ii1N
if (abs(H1(ii))>10e4) 防止H值太该频率信号放
XX1(ii)Yf(ii)*H1(ii)H(ii) 仿真结果
end
end
x1ifft(XX1)
subplot(121)
plot(ty1)
title(＇实际输出＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
subplot(122)
plot(treal(x1))
title(＇仿真输出＇)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
figure(10)
subplot(211)plot(ty2＇g＇)
xlabel(＇时间s＇)title(＇巴特沃思滤波器滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty1＇r＇)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
仿真短周期震仪短周期震仪窄带椭圆滤波器表示
ws[001 45]*2Fswp[01 38]*2Fs 通带阻带边界频率（化频率）
Rp1Rs20Nn512 通带波纹阻带衰减绘制频率特性数点数
[orderWn]ellipord(wpwsRpRs) 求取数字滤波器阶数化截止频率
[ba]ellip(orderRpRsWn) 阶数截止频率通带波纹阻带衰减设计滤波器
figure(11)
[H1f]freqz(baNnFs) 传递函数系数数点数采样频率求滤波器频率特性
subplot(211)plot(f20*log10(abs(H1)))
xlabel(＇频率Hz＇)ylabel(＇振幅dB＇)grid on
subplot(212)plot(f180pi*unwrap(angle(H1)))
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)grid on
figure(12)
y1filtfilt(baXt) 窄带滤波器输出
[H1f]freqz(baNFs＇whole＇) 震仪特性
XX1zeros(1N)
for ii1N 仿真结果
if (abs(H1(ii))>10e4)
XX1(ii)Yf(ii)*H1(ii)H(ii)
end
end
x1ifft(XX1)
plot(ty1treal(x1)＇r＇) 绘制输入信号
legend(＇实际输出＇＇仿真输出＇1)
xlabel(＇时间s＇)
ylabel(＇振幅＇)
grid on
仿真短周期震仪短周期震仪切雪夫2型滤波器表示
ws[001 45]*2Fswp[01 38]*2Fs
Rp1Rs20Nn512
[OrderWn]cheb2ord(wpwsRpRs)求取数字滤波器阶数化截止频率
[ba]cheby2(OrderRpWn)阶数截止频率通带波纹阻带衰减设计滤波器
figure(13)
[Hf]freqz(baNnFs)传递函数系数数点数采样频率求滤波器频率特性
y3filtfilt(baXt)
subplot(211)plot(f20*log10(abs(H)))画出宽带滤波器幅频特性
xlabel(＇＼lambda＇)ylabel(＇A(＼lambda)db＇)
title(＇宽频带滤波器幅频特性＇)
grid on
subplot(212)plot(fangle(H)) 画出宽带滤波器相频特性
xlabel(＇频率Hz＇)ylabel(＇相位^o＇)
title(＇宽频带滤波器相频特性＇)
grid on
figure(14)
subplot(211)plot(ty3＇g＇)
xlabel(＇时间s＇)title(＇切雪夫2型滤波器滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
subplot(212)plot(ty1＇r＇)
xlabel(＇时间s＇)title(＇椭圆宽带滤波器输出信号＇)
ylabel(＇振幅＇)
grid on
close allclear allclc
load hns1dat 读取数序列
Xthns1 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(1)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
title(＇P波＇)
xlabel(＇时间s＇)title(＇时间域＇)
title(＇P波＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf)) 绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
load hns2dat 读取数序列
Xthns2 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(2)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
title(＇S波＇)
xlabel(＇时间s＇)title(＇时间域＇)
title(＇S波＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf)) 绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
load hns3dat 读取数序列
Xthns3 数赋值变量
Fs50 设定采样率 单位(Hz)
dt1Fs 求采样间隔 单位(s)
Nlength(Xt) 序列长度
t[0N1]*dt 时间序列
Yffft(Xt) 信号进行快速Fourier变换(FFT)
figure(3)
subplot(211)plot([0N1]FsXt) 绘制原始值序列
title(＇面波＇)
xlabel(＇时间s＇)title(＇时间域＇)
title(＇面波＇)
grid on
subplot(212)plot([0N1]N*Fsabs(Yf))绘制信号振幅谱
xlabel(＇频率Hz＇)title(＇幅频图＇)
ylabel(＇振幅＇)
xlim([0 2]) 频率轴画出2Hz频率前部分
grid on
实验结果
图1 切雪夫1型宽频带滤波器椭圆宽带滤波器输出信号
面运动时域原始信号
图2 输入信号输出信号
图3 宽频带振幅相位
图4 短周期窄带振幅相位
图5 实际输出仿真输出
图6 巴特沃夫长周期实际输出仿真输出
图7 震波面波P波S波幅频图
图8 长周期窄带滤波器幅频特性
长周期窄带滤波器相频特性
六实验体会
通次实验进步复数字信号处理关滤波器基础解理实际身边处处数字信号处理相关知识应语音识采集处理震波观测直观证实数字信号处理门课程重性
次实验中实际操作中加强实践力巩固数字信号处理理知识培养解决实际问题力设计程中提高思考力动手力学理知识时明白应实际
次课程设计认识足认识学基础知识究竟运什领域运老师学耐心指导发现选择巴特沃斯切雪夫滤波器问题修改调试终应效果理实践相结合优势处受益匪浅
参考文献
[1]焦瑞莉 罗倩 汪毓铎 顾奕数字信号处理[M]机械工艺出版社pp184195
[2]httplgbougzcomcnhtmlautooldprotel99yyong1htm
袁宇波MatlabProtel设计微机保护中Butterworth模拟低通滤波器
[3]庆禹电力系统数字光电量测系统原理技术[J]电网技术200125(4) pp15
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容，不代表《香当网》观点或立场，请自行判断内容的真实性和可靠性！
该内容是文档的文本内容，更好的格式请下载文档
