1.5维普matlab程序

这是1.5维普的matlab程序

%--------------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