小海豚zc 发表于 2015-7-16 19:19

关于三分之一倍频程中有效值的求法问题

第一个问题:《MATLAB在振动信号处理中的应用》一书中有三分之一倍频程的程序(程序在下面),请高手帮我看看程序有没有问题,主要是中间求有效值那里
%三分之一倍频程处理
clear
%clc
%close all hidden
format long
%fni=input('三分之一倍频程处理-输入数据文件名','s');
%fid=fopen(fni,'r')
%sf=fscanf(fid,'%f',1); %读入采样频率值
%fno=fscanf(fid,'%d',1); %读入输出数据文件名
%x=fscanf(fid,'%f',); %读入输入数据存成列向量
%status=fclose(fid);
sf=200;
%fno='out6_4.mat';
%load y
%x=y;
N=20000;
nn=1:N;
t=nn./sf;
x=3*sin(2*pi*2.25*t)+3*sin(2*pi*15*t)+3*sin(2*pi*23*t);
%定义三分之一倍频程的中心频率
f=;
fc=;
%中心频率与下限频率的比值
oc6=2^(1/6);
nc=length(fc);
n=length(x);
nfft=2^nextpow2(n);
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
figure
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('有效值');
xlim()
grid on;
%fid=fopen(fno,'w');
%for k=1:m
% fprintf(fid,'%f%f\n',fc(k),yc(k));
%end
%status=fclose(fid);




xinglong-liu 发表于 2015-7-16 22:43

好像第41行“yc(j)=sqrt(var(real(b(1:n))));%计算每个中心频率段的有效值”改为“yc(j)=std(c)”即可

小海豚zc 发表于 2015-7-17 09:42

xinglong-liu 发表于 2015-7-16 22:43
好像第41行“yc(j)=sqrt(var(real(b(1:n))));%计算每个中心频率段的有效值”改为“yc(j)=std(c)”即可

首先谢谢你,还有几个问题想请教下
1.c是时域的,b是频域的,看本论坛的帖子说原理上都可以实现,那用b怎么做呢?
2.如果让输入是一个正弦信号,那均方根(有效值)应该是幅值除根2,可是用上面的程序却算不对(接近幅值除根2),误差是因为什么呢??是因为FFT的能量泄露吗??
3.有这样一种想法:既然信号可以做傅里叶级数展开,那我可以先把信号做FFT幅值谱,然后幅值谱除根2,再对每个频带做均方根。。这个做法可行吗??

TestGuru 发表于 2015-7-18 10:12

非整周期采样都会有泄漏问题,正弦信号用FFT做出来的谱线能量总会有泄漏,泄漏多少看所用的窗函数,能量总小于理论值。

小海豚zc 发表于 2015-7-18 10:59

TestGuru 发表于 2015-7-18 10:12
非整周期采样都会有泄漏问题,正弦信号用FFT做出来的谱线能量总会有泄漏,泄漏多少看所用的窗函数,能量总 ...

恩,有道理
能帮我看看楼上的第三个问题吗??用这种除根2求有效值的方式可行吗?

rossbin 发表于 2015-7-22 19:57

三分之一倍频程的求法在频域内是相应octave内的均方根值,频域内的计算是平方和开根号即可。

rossbin 发表于 2015-7-22 20:06

小海豚zc 发表于 2015-7-17 09:42
首先谢谢你,还有几个问题想请教下
1.c是时域的,b是频域的,看本论坛的帖子说原理上都可以实现,那用b ...

1、频域计算三分之一倍频程是相应octave内的rms值;
2、对简单的正弦信号,整周期采用,FFT到频谱后,频域的幅值需要scale,貌似要除以N,单边谱的话要除以N再乘以2。需要注意的是时域和频域的均方根值得算法。均方根值就是rms,时域先s,再m,最好r。但频域先s,然后直接r,不需要m,用matlab可以直接验证的。最近不做这个了,应当没有记错。
3、周期性信号才可展开傅里叶级数,非周期的要用傅里叶变换。个人认为不可行,你需要了解时域到频域信号的处理过程。

小海豚zc 发表于 2015-7-24 19:48

rossbin 发表于 2015-7-22 20:06
1、频域计算三分之一倍频程是相应octave内的rms值;
2、对简单的正弦信号,整周期采用,FFT到频谱后,频 ...
1.你说的octive是专门的一个软件吗?还是哪个概念?不好意思,我知道的太少了
2.你说的整周期采样,可以用上面的程序试验下,我自己试过了:一个正弦信号,即使整周期采样,还是达不到幅值除以根2(一个正选信号的有效值不就是幅值除以根2吗?),总是缺少一些,不清楚原因???还请麻烦您运行试下,帮我找下到底是什么问题

3.第三问题,其实我也一直有疑惑,有人说“功率信号才能做功率谱,能量信号才能做幅值谱”,可是功率谱有一种做法(Cooley-Turkey法)不就是在幅值谱的基础上调整幅值得到的吗?从程序算法的角度讲,这两种谱的求法一定存在关系,只不过从物理意义上讲就讲不通了。。   
然后回到我的问题,你说的非周期信号不能用fft,那非周期信号怎么用频域法?还有就是有限长度的非周期信号不是也可以拓展成周期信号吗?
不知道我说的对不对,还请您指正

rossbin 发表于 2015-7-24 19:54

小海豚zc 发表于 2015-7-24 19:48
1.你说的octive是专门的一个软件吗?还是哪个概念?不好意思,我知道的太少了
2.你说的整周期采样,可以 ...

你需要多看书。

小海豚zc 发表于 2015-7-24 22:56

rossbin 发表于 2015-7-24 19:54
你需要多看书。

哈哈,听君一席话,胜读十年书啊

我又想了下第三个问题:对一个信号,我没法保证对关注的频率都是整周期采样,所以,这样直接除根2 的算法确实是有问题的,

但是第二个问题还是考虑不清楚,原理上是对的

TestGuru 发表于 2015-7-25 17:08

整周期采样,算出来的正弦波的幅值一定会对的。否则检测程序是否有错误。

非整周期采样,需要采用适当的窗函数把往外泄漏的能量集中到主谱线附近,算出来的正弦波的幅度也一定会对的,否则检测程序是否有错误。

小海豚zc 发表于 2015-8-3 14:22

TestGuru 发表于 2015-7-25 17:08
整周期采样,算出来的正弦波的幅值一定会对的。否则检测程序是否有错误。

非整周期采样,需要采用适当的 ...

您帮忙看一下,我找不到问题

%三分之一倍频程处理
clear
%clc
%close all hidden
format long

sf=2000;
N=100000;
nn=1:N;
t=nn./sf;
x=2*sin(2*pi*2*t)+2*sqrt(2)*sin(2*pi*10*t)+2*sin(2*pi*30*t);
%定义三分之一倍频程的中心频率
f=;
fc=;
%中心频率与下限频率的比值
oc6=2^(1/6);
nc=length(fc);
n=length(x);
nfft=2^nextpow2(n);
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(c)));
end
figure
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('有效值');
xlim()
grid on;

GuoCL 发表于 2018-7-13 17:45

学习了,谢谢。

wwhblue 发表于 2020-2-17 17:24

fc最大8*10000,sf=200 ??? 不符合采用定理吧!?
页: [1]
查看完整版本: 关于三分之一倍频程中有效值的求法问题