我看了一些三点、五点求数值微分的程序，还是用符号计算的，这样用这些程序计算瞬时频率还是不靠谱吧，不 ...
忍不住想骂人！
为什么非得强调是符号计算的？强调符号计算时说明数值微分是连续的？那不用符号表示用啥表示计算过程？请楼主区分Matlab中符号计算和数值计算的意思。
不管你用什么差分算法算，都得不到和理论结果一致的结果。点数用多了也不见得会更好，相反如果变化大的情况下可能不如少数点。在时频工具箱中计算瞬时频率的程序默认是前后两点计算瞬时频率，当然也提供了参数选择多点计算瞬时频率，还可以做修改，比如引入AR过程对数据做预测做模型，但是不见得能得到很好的结果。
最后，不要混淆视听，实际过程中利用Matlab处理数据都是数值算法，都是近似的，只是近似程度高低的问题。 黄FAacos.m函数就是用acos函数求解的，但是里面特别考虑了反余弦后相位在0和pi附近的情况，程序里面是直接清除极大极小值之后再插值。所以直接acos离散获取相位会存在一些跳跃点的。 qqcd 发表于 2013-4-24 22:45 static/image/common/back.gif
黄FAacos.m函数就是用acos函数求解的，但是里面特别考虑了反余弦后相位在0和pi附近的情况，程序里面是直接清 ...
这样做是可以的，所以在利用diff做数值微分的时候，特别针对这种情况可以修改diff,添加代码部分处理当x=1时，确定其微分。而在楼主所提的问题中，符号计算不会出现x=1时其微分不等于0的情况。希望你将你提到的FAcos.m文件贴到后面！ %function = 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-
% = 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 normalizationfor an IMF ---loop A start
% 3.compute the phase angle by arc-cosine function
% 4.take difference for the phase angle
%2.Do the normalizationfor an IMF ---loop A end
%
function = FAacos(data, dt)
disp('WARNING: Applying the function to the residue can cause an error.');
%1.read the data,check input matrix
%----- Get dimensions
= size(data);
%----- Flip data if necessary
flipped=0;
if (npts < nimf)
flipped=1;
data=data';
= 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 normalizationfor 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 normalizationfor an IMF ---loop A end
end
%----- Flip again if data was flipped at the beginning
if (flipped)
f=f';
a=a';
end
qqcd 发表于 2013-4-24 22:45 static/image/common/back.gif
黄FAacos.m函数就是用acos函数求解的，但是里面特别考虑了反余弦后相位在0和pi附近的情况，程序里面是直接清 ...
这个函数我用过，效果还是不错的。不过我有个小疑问，这个函数找极大极小值后，得到的NaN个数与extr函数得到的个数不一样，以x=cos(16*pi*t+12*pi*t.*t)函数为例，黄程序里NaN个数是54个，而extr函数求得的x的极值点个数是27个，这是什么原因呢？ 学习了{:{51}:} {:{10}:}
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;
= hhspectrum(imf');
f=-f*fs;
figure(3)
plot(f(3,:))
xlabel('time');ylabel('frequence(Hz)');
~~~~~~~~~~~~
function = 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 =insfrequence(x,t)
sx=angle(x);
%计算瞬时频率
dt=diff(t);
dx=diff(sx);
f=(dx./dt)/(2*pi);
http://forum.vibunion.com/thread-126081-2-1.html
看这个帖子，这个问题在这个帖子里有解释。另外，发帖前请搜索主题。
页:
1
[2]