声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

12
返回列表 发新帖
楼主: shuihai707

[HHT] 符号和数值计算瞬时频率的差别

[复制链接]
发表于 2013-4-24 13:06 | 显示全部楼层

忍不住想骂人!
为什么非得强调是符号计算的?强调符号计算时说明数值微分是连续的?那不用符号表示用啥表示计算过程?请楼主区分Matlab中符号计算和数值计算的意思。
不管你用什么差分算法算,都得不到和理论结果一致的结果。点数用多了也不见得会更好,相反如果变化大的情况下可能不如少数点。在时频工具箱中计算瞬时频率的程序默认是前后两点计算瞬时频率,当然也提供了参数选择多点计算瞬时频率,还可以做修改,比如引入AR过程对数据做预测做模型,但是不见得能得到很好的结果。
最后,不要混淆视听,实际过程中利用Matlab处理数据都是数值算法,都是近似的,只是近似程度高低的问题。
回复 支持 反对
分享到:

使用道具 举报

发表于 2013-4-24 22:45 | 显示全部楼层
黄FAacos.m函数就是用acos函数求解的,但是里面特别考虑了反余弦后相位在0和pi附近的情况,程序里面是直接清除极大极小值之后再插值。所以直接acos离散获取相位会存在一些跳跃点的。
发表于 2013-4-25 16:29 | 显示全部楼层
qqcd 发表于 2013-4-24 22:45
黄FAacos.m函数就是用acos函数求解的,但是里面特别考虑了反余弦后相位在0和pi附近的情况,程序里面是直接清 ...

这样做是可以的,所以在利用diff做数值微分的时候,特别针对这种情况可以修改diff,添加代码部分处理当x=1时,确定其微分。而在楼主所提的问题中,符号计算不会出现x=1时其微分不等于0的情况。希望你将你提到的FAcos.m文件贴到后面!
发表于 2013-4-26 10:06 | 显示全部楼层
%function [f,a] = FAacos(data, dt)
%
% The function FAacos generates an arccosine frequency and amplitude
% of previously normalized data(n,k), where n is the length of the
% time series and k is the number of IMF components.
%
% FAacos finds the frequency by applying
% the arc-cosine function to the normalized data and
% checking the points where the slope of arc-cosine phase changes.
% Nature of the procedure suggests not to use the function
% to process the residue component.
%
% Calling sequence-
% [f,a] = FAacos(data, dt)
%
% Input-
%          data        - 2-D matrix of IMF components
%          dt            - time increment per point
% Output-
%          f            - 2-D matrix f(n,k) that specifies frequency
%          a            - 2-D matrix a(n,k) that specifies amplitude
%Note:
% Used by-
%        FA
%Reference:  
%
%written by
% Kenneth Arnold (NASA GSFC)        Summer 2003, Initial
%footnote: S.C.Su (2009/09/12)
%
%1.read the data,check input matrix
%2.Do the normalization  for an IMF ---loop A start
%   3.compute the phase angle by arc-cosine function   
%   4.take difference for the phase angle
%2.Do the normalization  for an IMF ---loop A end
%

function [f,a] = FAacos(data, dt)

disp('WARNING: Applying the function to the residue can cause an error.');

%1.read the data,check input matrix
    %----- Get dimensions
    [npts,nimf] = size(data);
   
    %----- Flip data if necessary
    flipped=0;
    if (npts < nimf)
        flipped=1;
        data=data';
        [npts,nimf] = size(data);
    end
   
    %----- Input is normalized, so assume that amplitude is always 1
    a = ones(npts,nimf);
   
    %----- Mark points that are above 1 as invalid (Not a Number)
    data(find(abs(data)>1)) = NaN;

%2.Do the normalization  for an IMF ---loop A start
    %----- Process each IMF
    for c=1:nimf
        %----- Compute "phase" by arc-cosine function
        acphase = acos(data(:,c));
        
%3.compute the phase angle by arc-cosine function
    %----- Mark points where slope of arccosine phase changes as invalid
    %after difference the arc-cosine function will be discontinuous at those points
      for i=2:npts-1
          prev = data(i-1,c);%value of 'i-1' position
          cur = data(i,c);   %value of  'i'  position
          next = data(i+1,c);%value of 'i+1' position
         
          %Check for local max and local min, set them as NaN
          if (prev < cur & cur > next) | (prev > cur & cur < next)
              acphase(i) = NaN;
              
          end
      end
           
      clear prev cur next
       %clear syntax is important for this kind of algorithm-S.C.Su
       %choose some value with out clear makes chaos
%4.take difference for the phase angle                                 
    %----- Get phase differential frequency--calculate I.F. values
    acfreq = abs(diff(acphase))/(2*pi*dt);
   
    %----- Mark points with negative frequency as invalid
    acfreq(find(acfreq<0)) = NaN;
   
    %----- Fill in invalid points using a spline
    legit = find(~isnan(acfreq));%legit=allvalues-NaN
      %interpolate for NaN positions  
      if (length(legit) < npts)
        f(:,c) = spline(legit, acfreq(legit), 1:npts)';
      else
        f(:,c) = acfreq;
      end
%2.Do the normalization  for an IMF ---loop A end
    end

%----- Flip again if data was flipped at the beginning
if (flipped)
    f=f';
    a=a';
end

评分

1

查看全部评分

 楼主| 发表于 2013-4-26 12:48 | 显示全部楼层
qqcd 发表于 2013-4-24 22:45
黄FAacos.m函数就是用acos函数求解的,但是里面特别考虑了反余弦后相位在0和pi附近的情况,程序里面是直接清 ...

这个函数我用过,效果还是不错的。不过我有个小疑问,这个函数找极大极小值后,得到的NaN个数与extr函数得到的个数不一样,以x=cos(16*pi*t+12*pi*t.*t)函数为例,黄程序里NaN个数是54个,而extr函数求得的x的极值点个数是27个,这是什么原因呢?
发表于 2013-4-26 15:14 | 显示全部楼层
学习了
发表于 2013-4-27 16:58 | 显示全部楼层
发表于 2013-10-10 17:43 | 显示全部楼层

HHT处理一般信号,为什么会出现如此问题

本人以前用网上提供的HHT程序来分析桥梁结构在移动荷载下的某点响应,程序老出错,比如infreq程序得到的结果错误,所以自己在人家的基础上大改了一下。现在就用改好的先试试简单的几个信号。分析x1=sin(t)和x2=sin(2*pi*t的叠加信号,得到如图结果,为什么瞬时频率图怎么是这样的呢?应该是直线才对啊,为什么会周期性的出现极值呢?请指点,能否提供一整套HHT的程序?谢谢nzhenhua83@gmail.com
图和程序如下
t=1:0.01:100;
fs=100;
x1=sin(t);                  %f=1/(2*pi)=0.1629Hz
x2=sin(2*pi*t);          %f=1Hz
Y=x1+x2;

Nstd=0;NE=1;
allmode=eemd(Y,Nstd,NE);   

imf=allmode;
[A,f,tt] = hhspectrum(imf');

f=-f*fs;
figure(3)
plot(f(3,:))
xlabel('time');ylabel('frequence(Hz)');

~~~~~~~~~~~~
function [A,f,tt] = hhspectrum(imf,t,l)
if nargin < 2
t=1:size(imf,2);

end
if nargin < 3
  l=1;
end
if nargin < 4
  aff = 0;
end
for i=1:(size(imf,1)-1)
  an(i,:)=hilbert(imf(i,:)')';
  f(i,:)=insfrequence(an(i,:),t)';
  A=abs(an(:,l+1:end-l));
end

~~~~~~~~~~~
function [f]=insfrequence(x,t)
sx=angle(x);
%计算瞬时频率
dt=diff(t);
dx=diff(sx);
f=(dx./dt)/(2*pi);

imf3瞬时频率

imf3瞬时频率

imf2瞬时频率

imf2瞬时频率

EMD分解信号

EMD分解信号

原始信号

原始信号
发表于 2013-10-10 20:04 | 显示全部楼层
http://forum.vibunion.com/thread-126081-2-1.html
看这个帖子,这个问题在这个帖子里有解释。另外,发帖前请搜索主题。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 05:28 , Processed in 0.064519 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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