Halcom 发表于 2017-8-5 10:50:59

来自三个方向的信号强度计算

SignalEW、SignalNS、SignalUD不同方向的信号,融合求解强度Intensity
软件:MATLAB
function [ IntensityV ] = IntensityCalFunc(SignalEW, SignalNS, SignalUD)

MaxEW=max(abs(SignalEW));
MaxNS=max(abs(SignalNS));
MaxUD=max(abs(SignalUD));
SignalLength=length(SignalEW);

%FFT
Fs = 50;                  % Sampling frequency 50Hz
T = 1./Fs;                  % Sample time
L = length(SignalEW);       % Length of signal
t = (0:L-1)*T;            % Time vector

%FFT Start
SignalEWfft = fft(SignalEW);
SignalNSfft = fft(SignalNS);
SignalUDfft = fft(SignalUD);

% fft scale
SampleSpace=1./Fs;
f_val=1./(L*SampleSpace);
f_results=zeros(L,1);
f_results=f_results';

f_N=floor(L/2);
f_p1=0:(f_N-1);
f_results(1:f_N)=f_p1;
f_p2=(-floor(L/2)):-1;
f_results(f_N+1:end)=f_p2;
f=abs(f_results.*f_val);

%make Filter
% period filter window
winX = sqrt(1./f(2:end));
winX=horzcat(0,winX);
%high cut filter
y = f / 10.;
winY = 1./sqrt(1.0 + 0.694*power(y,2) + 0.241*power(y,4) + 0.0557*power(y,6) + 0.009664*power(y,8) + 0.00134*power(y,10) + 0.000155*power(y,12));
%low cut filter
winZ = sqrt(1.0 - exp( power( -1.0 * (f / 0.5), 3) ));
%three windows composition
win1=winX.*winY;
win = win1.*winZ;

% signal fft filtering
spectrum_winE = win.*SignalEWfft ;
spectrum_winN = win.*SignalNSfft;
spectrum_winU = win.*SignalUDfft;

%fourier inverse
resyn_sigE = ifft(spectrum_winE);
resyn_sigN = ifft(spectrum_winN);
resyn_sigU = ifft(spectrum_winU);

%wave composition
S2 = sqrt(power(resyn_sigE,2) + power(resyn_sigN,2) + power(resyn_sigU,2));

%over a (0.3sec)
maxS2 = sort(S2, 'descend');
I = 2*log10(maxS2(int32(0.3*Fs))) + 0.94;
real(I);
IntensityV= ( round(real(I),2) * 10 ) / 10.0;
end
参考:
【1】http://halcom.cn/forum.php?mod=viewthread&tid=2928&extra=page%3D1
【2】一维信号小波去噪,傅里叶幅值和相位角计算
【3】去随机飘逸噪声



页: [1]
查看完整版本: 来自三个方向的信号强度计算