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

Halcom

SignalEW、SignalNS、SignalUD不同方向的信号，融合求解强度Intensity
软件：MATLAB

1. function [ IntensityV ] = IntensityCalFunc(SignalEW, SignalNS, SignalUD)
   
2. 
   
3. MaxEW=max(abs(SignalEW));
   
4. MaxNS=max(abs(SignalNS));
   
5. MaxUD=max(abs(SignalUD));
   
6. SignalLength=length(SignalEW);
   
7. 
   
8. %FFT
   
9. Fs = 50;                    % Sampling frequency 50Hz
   
10. T = 1./Fs;                  % Sample time
    
11. L = length(SignalEW);       % Length of signal
    
12. t = (0:L-1)*T;              % Time vector
    
13. 
    
14. %FFT Start
    
15. SignalEWfft = fft(SignalEW);
    
16. SignalNSfft = fft(SignalNS);
    
17. SignalUDfft = fft(SignalUD);
    
18. 
    
19. % fft scale
    
20. SampleSpace=1./Fs;
    
21. f_val=1./(L*SampleSpace);
    
22. f_results=zeros(L,1);
    
23. f_results=f_results';
    
24. 
    
25. f_N=floor(L/2);
    
26. f_p1=0:(f_N-1);
    
27. f_results(1:f_N)=f_p1;
    
28. f_p2=(-floor(L/2)):-1;
    
29. f_results(f_N+1:end)=f_p2;
    
30. f=abs(f_results.*f_val);
    
31. 
    
32. %make Filter
    
33. % period filter window
    
34. winX = sqrt(1./f(2:end));
    
35. winX=horzcat(0,winX);
    
36. %high cut filter
    
37. y = f / 10.;
    
38. 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));
    
39. %low cut filter
    
40. winZ = sqrt(1.0 - exp( power( -1.0 * (f / 0.5), 3) ));
    
41. %three windows composition
    
42. win1=winX.*winY;
    
43. win = win1.*winZ;
    
44. 
    
45. % signal fft filtering
    
46. spectrum_winE = win.*SignalEWfft ;
    
47. spectrum_winN = win.*SignalNSfft;
    
48. spectrum_winU = win.*SignalUDfft;
    
49. 
    
50. %fourier inverse
    
51. resyn_sigE = ifft(spectrum_winE);
    
52. resyn_sigN = ifft(spectrum_winN);
    
53. resyn_sigU = ifft(spectrum_winU);
    
54. 
    
55. %wave composition
    
56. S2 = sqrt(power(resyn_sigE,2) + power(resyn_sigN,2) + power(resyn_sigU,2));
    
57. 
    
58. %over a (0.3sec)
    
59. maxS2 = sort(S2, 'descend');
    
60. I = 2*log10(maxS2(int32(0.3*Fs))) + 0.94;
    
61. real(I);
    
62. IntensityV= ( round(real(I),2) * 10 ) / 10.0;
    
63. end
    

复制代码

参考：
【1】http://halcom.cn/forum.php?mod=viewthread&tid=2928&extra=page%3D1
【2】一维信号小波去噪，傅里叶幅值和相位角计算
【3】去随机飘逸噪声