matlab编程遇到的问题

海峡的风

运行之后发现fun为一实数，不是含有未知数k的式子，所以它的一阶导数fk1,二阶导数fk2均是0，请各位高手指点一下，问题出在哪，该怎么改呢（真的很急啊，谢谢大家帮忙）
clc
clear
syms k;
N=1024*4;
Fs=2000;
t=(0:N-1)/Fs;
x=10*(1+sin(2*pi*10*t)).*sin(2*pi*200*t+8*sin(2*pi*10*t))+10*sin(2*pi*400*t);
[c]=mystfta(x,32,16,Fs);%短时傅立叶变化
ac=abs(c);
[u,v]=size(ac);
SK=[];
  for kk=1:u
    y=ac(kk,:);
    Temp=kurtosis(y);%计算谱峭度
    SK=[SK Temp];
    x1=[SK];
  end
f1=(0:u-1)/u*Fs/2;
figure;
subplot(311)
plot(x);           %绘制原信号
subplot(312)
plot(f1,x1,'r');  %做出鞘度谱图
%滤波设计
n=[SK];
for i=1:u
   a=n(1,i);
   M=sqrt(a/k);  %维纳滤波器的频率响应函数
   y=x.*M;         %输出信号
   fun=sum(y.*y.*y.*y)/(sum(y.*y)*sum(y.*y))-2;  %输出信号的峭度
   %fun=kurtosis(y);
   fk1=diff(fun,'k');   %使用牛顿法求最优解
   fk2=diff(fk1,'k');
   k=a;
for n=1:100
    f0=subs(fun);
    ff1=subs(fk1);  
    ff2=subs(fk2);
   if abs(ff1)<= 0.001
     
      k=vpa(k,4);
      f0=vpa(f0,4);
      kkk(i)=k; %使得导数接近于0的k值
      ff0(i)=f0;
      [m,n]
   else
      K=k-ff1/ff2;
      k=K;
      
   end
end
end
[m,n]=min(ff0);
k=K(n);%得到最优解
有一个函数fmincon 可以用来求最优解 ，我也试了一下，但也有问题说这个函数用的不对，大伙能找出原因么 ，先谢谢大家了
clc
clear
N=1024*4;
Fs=2000;
t=(0:N-1)/Fs;
x=10*(1+sin(2*pi*10*t)).*sin(2*pi*200*t+8*sin(2*pi*10*t))+10*sin(2*pi*400*t);
    [c]=mystfta(x,32,16,Fs);
    ac=abs(c);
    [u,v]=size(ac);
    SK=[];
  for kk=1:u
    y=ac(kk,:);
    %Temp=sum(y.*y.*y.*y)/(sum(y.*y)*sum(y.*y))-3;
    Temp=kurtosis(y);
    SK=[SK Temp];
    x1=[SK];
  end
f1=(0:u-1)/u*Fs/2;
figure;
%f=(0:)/Nw*Fs/2
subplot(311)
plot(x);
subplot(312)
plot(f1,x1,'r');
%维纳滤波设计
n=[SK];
for i=1:u
    syms k;
        a=n(1,i);
        M=sqrt(a/k);
        y=x.*M;
        %fun='sum(y.*y.*y.*y)/(sum(y.*y)*sum(y.*y))-3';
        [k,fval]=fmincon(@(k) sum(y.*y.*y.*y)/(sum(y.*y)*sum(y.*y))-3,a,-1,-a,[],[],a,inf);
      %[K,fval]= fseminf(fun,a,1,@mycon);
        K(i)=k;
        Wmin(i)=fval;
      
end
[a,b]=min(Wmin);
k=K(b);
y=x*sqrt(SK./k);
subplot(313)
plot(y);
其中mystfta.m文件为：
function [c] =mystfta(x,shift,wlen,fs)
% x--signal to be processing
% shift--shift length of moving window
% wlen--window length
% window--type of window
s=size(x);
if s(2)>1
   x=x';
end
lv=max(size(x));  % equ. length(x)
x=x-mean(x);
n=(lv-wlen)/shift+1;

w=hanning(wlen);
deltf=fs/wlen;
deltt=1./fs;
c = [];
for k = 1:n
        i=(k-1)*shift+1;
        s=x(i:i+wlen-1);
        s=s.*w;
        s = fft(s);
        m=length(s);
        mm=fix(m/2);
        c = [c s(1:mm,:)];
end
rc=real(c);
ic=imag(c);
%ac=abs(c(1:wlen/2,:));
save Rcfd rc -ascii;
save Icfd ic -ascii;

hehehejunfei

您好，请问这个问题是怎么解决的？

xjchina

绑定绑定绑定绑定matlab在振动信号处理中的应用