声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1974|回复: 2

[编程技巧] matlab编程遇到的问题

[复制链接]
发表于 2009-8-30 21:49 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
运行之后发现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;
回复
分享到:

使用道具 举报

发表于 2014-8-23 11:05 | 显示全部楼层
您好,请问这个问题是怎么解决的?
发表于 2014-8-25 20:37 | 显示全部楼层
绑定绑定绑定绑定matlab在振动信号处理中的应用
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-13 10:42 , Processed in 0.075506 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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