三分之一倍频程程序
更新时间:2023-11-01 05:18:01 阅读量: 综合文库 文档下载
- 三分野推荐度:
- 相关推荐
方法一:%A计权声压级频谱分析 clc; clear; close all;
y=wavread('abc.wav');
fs=51200;%采样频率 p0=2e-5;%参考声压
f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率 f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%% -16000Hz A声级计权值
cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6]; t1=1; t2=2;
x=y(t1*fs+1:t2*fs);%截取需要处理的数据段 n=length(x);
t=(0:1/fs:(n-1)/fs);
subplot(221);
plot(t,x);%瞬时声压时程图
w=hanning(n); %汉宁窗
xx=1.633*x.*w; %加汉宁窗(恢复系数为1.633)
nfft=2^nextpow2(n);
%nextpow2(n)-取最接近的较大2次幂 a = fft(xx,nfft);
f = fs/2*linspace(0,1,nfft/2);
w=2*abs(a(1:nfft/2)/n); subplot(222);
plot(f,w);%绘制频谱图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %1/3倍频程计算 oc6=2^(1/6); nc=length(cf);
%下面这个求1/3倍频程的程序是按照振动振级计算那个来的 for j=1:nc
fl=fc(j)/oc6; fu=fc(j)*oc6;
nl=round(fl*nfft/fs+1); nu=round(fu*nfft/fs+1); if fu>fs/2 m=j-1; break; end
b=zeros(1,nfft); b(nl:nu)=a(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1); c=ifft(b,nfft);
yc(j)=sqrt(var(real(c(1:n)))); end
aj_sumn=0; for i=1:nc
Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级 end
%%%%% for jj=1:nc
aj_sumn=aj_sumn+10^(0.1*Lp1(j)); end
Lp=10*log10(aj_sumn);%未计权总声压级
subplot(223);%绘制未计权1/3倍频程声压级图谱
bar(Lp1(1:nc)); gg=zeros(1,nc); for i=1:nc
gg(1:nc)=fc(1:nc); end
ggg=1:nc; set(gca,'xtick',ggg); set(gca,'xticklabel',gg);
%%%%%A计权1/3倍频程声压级 Lap=Lp1+cf; aj_sum=0; for j=1:nc
aj_sum=aj_sum+10^(0.1*Lap(j)); end
LA=10*log10(aj_sum);a计权总声压级
subplot(224);%绘制A计权1/3倍频程声压级图谱 bar(Lap(1:nc)); gg=zeros(1,nc); for i=1:nc
gg(1:nc)=fc(1:nc); end
ggg=1:nc; set(gca,'xtick',ggg); set(gca,'xticklabel',gg); 方法二: clc; clear; close all;
%时域分析
y=wavread('abc.wav'); %频域分析
fs=51200;%采样频率 p0=2e-5;%参考声压
f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率 f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%% -16000Hz A声级计权值
cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6];
n=length(y);
t=(0:1/fs:(n-1)/fs); h1=figure; plot(t,y);
title('瞬时声压时程'); xlabel('Time(s)');
ylabel('Sound Presure Value(Pa)'); % t1=0; t2=4;
x=y(t1*fs+1:t2*fs);%截取需要处理的数据段 n=length(x);
%t=(0:1/fs:(n-1)/fs);
%plot(t,x);%瞬时声压时程图
w=1.633*hanning(n); %汉宁窗(恢复系数为1.633)
%w=1.812*blackmanharris(n); %布拉克曼窗(功率相等恢复系数1.812) xx=x.*w; %加汉宁窗
nfft=2^nextpow2(n); %nextpow2(n)-取最接近的较大2次幂 a = fft(xx,nfft)/n;
%f = fs/2*linspace(0,1,nfft/2); w=2*abs(a(1:nfft/2));
oc6=2^(1/6); nc=length(cf);
for j=1:nc
fl=fc(j)/oc6; fu=fc(j)*oc6;
nl=round(fl*nfft/fs+1); nu=round(fu*nfft/fs+1); if fu>fs/2 m=j-1; break; end
p=w(nl:nu); lp=length(p); k=0;
for ii=1:lp if ii+2>lp break end
if p(ii+1)>p(ii)&&p(ii+1)>p(ii+2) k=k+1;
pp(k)=p(ii+1)/sqrt(2); end end
p2(j)=sum(pp.*pp); end
Lp=10*log10(p2/p0^2); for jj=1:length(Lp)
Lp1(jj)=10^(Lp(jj)/10); end
Lpt=10*log10(sum(Lp1))
h2=figure;
mm=nc;
bar(Lp(1:mm)); gg=zeros(1,mm); for i=1:mm
gg(1:mm)=fc(1:mm); end
ggg=1:mm; set(gca,'xtick',ggg); set(gca,'xticklabel',gg);
set(gcf,'PaperPosition',[1,1,40,20]) set(gca,'fontsize',10) xlabel('Frequency(Hz)'); ylabel('SPL(dB)'); title('未计权声压级'); grid on;
方法三:
% 三分之一倍频程处理 clear; clc;
close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s = xlsread('ay.xls');%输入时程数据 sf=256; %采样频率
x=s(:,2); %定义三分之一倍频程的中心频率
f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00]; fc=[f,10*f,100*f,1000*f,10000*f]; %中心频率与下限频率的比值 oc6=2^(1/6);
%取中心频率总的长度 nc=length(fc);
%输入数据的长度 n=length(x);
%大于并接近n的2的幂次方长度 nfft=2^nextpow2(n); ?T变换 a=fft(x,nfft); for j=1:nc %下线频率 fl=fc(j)/oc6; %上限频率 fu=fc(j)*oc6;
%下限频率对应的序号 nl=round(fl*nfft/sf+1);
%上限频率对应的序号 nu=round(fu*nfft/sf+1);
%如果上相频率大于折叠频率则循环中断 if fu>sf/2 m=j-1;break end
%以每个中心频率段为通带进行带通频率滤波 b=zeros(1,nfft); b(nl:nu)=a(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1); c=ifft(b,nfft);
%计算对应每个中心频段的有效值 yc(j)=sqrt(var(real(b(1:n)))); end
%绘制输入时程曲线图形 subplot(2,1,1); t=0:1/sf:(n-1)/sf; plot(t,x);
xlabel('时间(s)'); ylabel('加速度(g)'); grid on;
%绘制三分之一倍频程有效值图形 subplot(2,1,2);
plot(fc(1:m),yc(1:m));
xlabel('频率(Hz)'); ylabel('有效值'); grid on;
%保存倍频程数据 fid=fopen(fno,’w’); for k=1:m;
fprintf(fid,’%f %f\\n’,fc(k),yc(k)); end
status=fclose(fid);
%上限频率对应的序号 nu=round(fu*nfft/sf+1);
%如果上相频率大于折叠频率则循环中断 if fu>sf/2 m=j-1;break end
%以每个中心频率段为通带进行带通频率滤波 b=zeros(1,nfft); b(nl:nu)=a(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1); c=ifft(b,nfft);
%计算对应每个中心频段的有效值 yc(j)=sqrt(var(real(b(1:n)))); end
%绘制输入时程曲线图形 subplot(2,1,1); t=0:1/sf:(n-1)/sf; plot(t,x);
xlabel('时间(s)'); ylabel('加速度(g)'); grid on;
%绘制三分之一倍频程有效值图形 subplot(2,1,2);
plot(fc(1:m),yc(1:m));
xlabel('频率(Hz)'); ylabel('有效值'); grid on;
%保存倍频程数据 fid=fopen(fno,’w’); for k=1:m;
fprintf(fid,’%f %f\\n’,fc(k),yc(k)); end
status=fclose(fid);
正在阅读:
三分之一倍频程程序11-01
中国人民解放军各集团军编制战斗序列大全05-02
危险源辨识与风险评价一览表01-05
五年级分段收费问题10-24
《建筑制图》练习题03-17
特种设备管理人员复习题(A3) —压力管道安全管理18710-07
五行属土的字姓名学解释06-01
传媒大学外国文学史真题考点分布修订稿05-03
工作中激励的经典语录02-07
- 多层物业服务方案
- (审判实务)习惯法与少数民族地区民间纠纷解决问题(孙 潋)
- 人教版新课标六年级下册语文全册教案
- 词语打卡
- photoshop实习报告
- 钢结构设计原理综合测试2
- 2014年期末练习题
- 高中数学中的逆向思维解题方法探讨
- 名师原创 全国通用2014-2015学年高二寒假作业 政治(一)Word版
- 北航《建筑结构检测鉴定与加固》在线作业三
- XX县卫生监督所工程建设项目可行性研究报告
- 小学四年级观察作文经典评语
- 浅谈110KV变电站电气一次设计-程泉焱(1)
- 安全员考试题库
- 国家电网公司变电运维管理规定(试行)
- 义务教育课程标准稿征求意见提纲
- 教学秘书面试技巧
- 钢结构工程施工组织设计
- 水利工程概论论文
- 09届九年级数学第四次模拟试卷
- 三分
- 之一
- 程序
- 新建 每次经历都是一次成长
- 有机化学习题答案
- 大学物理II练习册答案7
- 2019年高考主观创新题得体 答案
- 民俗文化的旅游开发
- (正式)淡发展中的中国股市过度投机性
- 2013 年深圳市高三年级第二次调研考试理科综合 - 图文
- 岗位安全标准化作业指导书模板
- 申论55需注意的
- 当代画竹子最有名画家 花鸟画竹子作品欣赏 - 图文
- 幼儿园游戏
- 浙江省稽阳市联谊学校2016届高三下学期4月联考地理试卷(解析) - 图文
- 高考各科成绩至少提高10分的小细节
- 江苏省高层次创新创业人才引进计划申报书(A类)
- 美岱召镇境内旅游资源分布情况说明
- 第4章需求开发与需求管理
- RADVISION深度测试流程
- 糖类 - 图文
- 2011-2012江苏省如东县七年级下英语试卷
- 前男友另觅新欢还来联系你,这样的情况该怎么应对?