频谱瀑布图_matlab例程_

clear; clc; close all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %以下为输入% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t1=0;%输入起始时间，若为0表示默认从0开始，若t1输入1表示从1s开始 t2=10;%输入终止时间，若为0表示分析全部，若t2输入2表示到2s结束 str='C:\Users\Administrator\Desktop\关键件 大族\大族j2 0629.wav';%文件路径和文件名 range1=16;%%绘图下截止频率贴近 2 range=5000;%%绘图上截止频率 inf=500;%%fft分析上限 inf2=50;%%软件读取上限 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %以上为输入% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [x1,Fs]=audioread(str); if t1==0 t1=t1+1; else t1=t1*Fs; end if t2==0 t2=length(x1); else t2=t2*Fs; end samples=[t1,t2]; [x1,Fs] = audioread(str,samples); x=x1(:,1)'; figure(1) subplot(211); t=(0:length(x)-1)/Fs; plot(t,x); title('时域分析图'); subplot(212); y=fft(x,Fs); df=Fs/(length(y)-1); fx=df*(0:length(y)-1); plot(fx,abs(y)); hold on plot([0,range],[200,200],'r--','LineWidth',2); % 绘制y=200的直线 plot([0,range],[150,150],'g--','LineWidth',2); % 绘制y=100的直线 axis([range1 range 0 inf]); title('频域分析图'); Pref=2e-5; % 参考声压 % 1/3倍频程滤波器中心频率 ff = [ 20, 25 31.5 40, 50 63 80, 100 125 160,... 200 250 315, 400 500 630, 800 1000 1250, 1600 2000 2500, ... 3150 4000 5000, 6300 8000 10000, 12500 16000]; nc=length(ff); % 1/3滤波器个数 P = zeros(1,nc); % 初始化 m = length(x); % x的长度 oc6=2^(1/6); % 倍频程的比例 for i=1:nc fl=ff(i)/oc6; % 求出1/3倍频程低截止频率 fu=ff(i)*oc6; % 求出1/3倍频程高截止频率 % 调用fdesign+designign函数计算滤波器系数集合Hd d=fdesign.bandpass('N,F3DB1,F3DB2',8,fl,fu,Fs); Hd=design(d); y = filter(Hd,x); % 滤波 P(i) = sum(y.^2)/m; % 计算输出信号的均方值 end % 计算各频带的声压级和总声压级 Psum=0; for i=1 : nc Pow(i) = 10*log10(P(i)/Pref^2);% 计算各频带的声压级 Psum=Psum+P(i); % 能量相加 end Lps=10*log10(Psum/Pref^2); % 计算总声压级 fprintf('LPS=%5.6fdB\n',Lps); figure(2) bar(Pow); set(gca,'XTick',[2:3:30]); grid set(gca,'XTickLabels',ff(2:3:length(ff))); xlabel('中心频率/Hz'); ylabel('声压级/dB'); title('三分之一倍频程滤波器输出声压级') set(gcf,'color','w'); figure(3) specgram(x1,2048,Fs,2048,1536);