请选择 进入手机版 | 继续访问电脑版

Hello Mat

 找回密码
 立即注册
查看: 4055|回复: 0

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

[复制链接]

1064

主题

1239

帖子

13

金钱

管理员

Rank: 9Rank: 9Rank: 9

积分
21393
发表于 2017-8-5 10:50:59 | 显示全部楼层 |阅读模式
SignalEW、SignalNS、SignalUD不同方向的信号,融合求解强度Intensity
软件:MATLAB
  1. function [ IntensityV ] = IntensityCalFunc(SignalEW, SignalNS, SignalUD)

  2. MaxEW=max(abs(SignalEW));
  3. MaxNS=max(abs(SignalNS));
  4. MaxUD=max(abs(SignalUD));
  5. SignalLength=length(SignalEW);

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

  11. %FFT Start
  12. SignalEWfft = fft(SignalEW);
  13. SignalNSfft = fft(SignalNS);
  14. SignalUDfft = fft(SignalUD);

  15. % fft scale
  16. SampleSpace=1./Fs;
  17. f_val=1./(L*SampleSpace);
  18. f_results=zeros(L,1);
  19. f_results=f_results';

  20. f_N=floor(L/2);
  21. f_p1=0:(f_N-1);
  22. f_results(1:f_N)=f_p1;
  23. f_p2=(-floor(L/2)):-1;
  24. f_results(f_N+1:end)=f_p2;
  25. f=abs(f_results.*f_val);

  26. %make Filter
  27. % period filter window
  28. winX = sqrt(1./f(2:end));
  29. winX=horzcat(0,winX);
  30. %high cut filter
  31. y = f / 10.;
  32. 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));
  33. %low cut filter
  34. winZ = sqrt(1.0 - exp( power( -1.0 * (f / 0.5), 3) ));
  35. %three windows composition
  36. win1=winX.*winY;
  37. win = win1.*winZ;

  38. % signal fft filtering
  39. spectrum_winE = win.*SignalEWfft ;
  40. spectrum_winN = win.*SignalNSfft;
  41. spectrum_winU = win.*SignalUDfft;

  42. %fourier inverse
  43. resyn_sigE = ifft(spectrum_winE);
  44. resyn_sigN = ifft(spectrum_winN);
  45. resyn_sigU = ifft(spectrum_winU);

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

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



算法QQ  3283892722
群智能算法链接http://halcom.cn/forum.php?mod=forumdisplay&fid=73
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Python|Opencv|MATLAB|Halcom.cn

GMT+8, 2022-7-3 13:53 , Processed in 0.174347 second(s), 21 queries .

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表