%--------------1又1/2維譜（雙譜在平面f1=f2上的投影）-----------------
clear?all;close?all;clc;
%time=xlsread（‘C:\Documents?and?Settings\Administrator\桌面\匯總插值數據.xls‘‘震蕩圖1‘‘A1:A2001‘）;%取時間值
CurrentVal=xlsread（‘C:\Documents?and?Settings\Administrator\桌面\匯總插值數據.xls‘‘故障韓‘‘B1:B2001‘）;%取電流值
subplot（211）;plot（CurrentVal）;title（‘波形一原始波形‘）;xlabel（‘時間/ms‘）;ylabel（‘電流幅值/V‘）;grid;
N=length（CurrentVal）;
%k=20;%把數據分成4段
%m=100;%每段數據500個點
%樣點數
%N=128*64;
%采樣頻率
fs=10000;
%頻率參數
%f1=0.6381;
%f2=0.8345;
%f3=f1+f2;
%f4=0.4909;
%f5=1.7671;
%f6=2;
%仿真信號
%a=2*pi*rand（5）;
%x=cos（2*pi*f1/fs*（1:N）+a（1））+cos（2*pi*f2/fs*（1:N）+a（2））+...
??%cos（2*pi*f3/fs*（1:N）+a（1）+a（2））+cos（2*pi*f4/fs*（1:N）+a（4））+...
??%cos（2*pi*f5/fs*（1:N）+a（5））+cos（2*pi*f6/fs*（1:N）+a（4）+a（5））;
%估計三階累積量
K=20;??%分成64段
M=fix（N/K）;?%每段的樣點數
%分段，去均值
xk=zeros（KM）;
for?i=0:K-1
????y=x（i*M+1:（i+1）*M）;
????y=y-mean（y）;
????xk（i+1:）=y;
end
nfft=256;
%計算累積量
c=zeros（Knfft）;
for?i=0:K-1
??for?t=-（nfft/2-1）:nfft/2
????s1=max（11-t）;%%這個地方不知道怎么修改
????s2=min（MM-t）;%%
????for?n=s1:s2
????j=xk（i+1n）*xk（i+1n+t）*xk（i+1n+t）;
????c（i+1t+nfft/2）=sum（j）*（1/M）;%%必要的時候自己用簡單的序列來測試下
????end
??end
??
end
cc=（1/K）*sum（c）;
%FFT變換
Y=fft（ccnfft）;
P=Y.*conj（Y）/nfft;
f=fs*（0:nfft/2-1）/nfft;
plot（fP（1:nfft/2））;
xlabel（‘頻率‘）;
ylabel（‘幅值‘）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????1463??2013-01-15?13:44??1.5.m

-----------?---------??----------?-----??----

?????????????????1463????????????????????1