怎样在MATLAB中划出一个函数的包络线

Posted 2023-05-08

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了怎样在MATLAB中划出一个函数的包络线相关的知识，希望对你有一定的参考价值。

下面是一系列关于MATLAB的包络线的程序：
%这是定义了一个函数：
function [up,down] = envelope(x,y,interpMethod)

%ENVELOPE gets the data of upper and down envelope of the known input (x,y).
%
% Input parameters:
% x the abscissa of the given data
% y the ordinate of the given data
% interpMethod the interpolation method
%
% Output parameters:
% up the upper envelope, which has the same length as x.
% down the down envelope, which has the same length as x.
%
% See also DIFF INTERP1

% Designed by: Lei Wang, , 11-Mar-2003.
% Last Revision: 21-Mar-2003.
% Dept. Mechanical & Aerospace Engineering, NC State University.
% $Revision: 1.1 $ $Date: 3/21/2003 10:33 AM $

if length(x) ~= length(y)
error('Two input data should have the same length.');
end

if (nargin < 2)|(nargin > 3),
error('Please see help for INPUT DATA.');
elseif (nargin == 2)
interpMethod = 'linear';
end

% Find the extreme maxim values
% and the corresponding indexes
%----------------------------------------------------
extrMaxValue = y(find(diff(sign(diff(y)))==-2)+1);
extrMaxIndex = find(diff(sign(diff(y)))==-2)+1;

% Find the extreme minim values
% and the corresponding indexes
%----------------------------------------------------
extrMinValue = y(find(diff(sign(diff(y)))==+2)+1);
extrMinIndex = find(diff(sign(diff(y)))==+2)+1;

up = extrMaxValue;
up_x = x(extrMaxIndex);

down = extrMinValue;
down_x = x(extrMinIndex);

% Interpolation of the upper/down envelope data
%----------------------------------------------------
up = interp1(up_x,up,x,interpMethod);
down = interp1(down_x,down,x,interpMethod); 参考技术A 下面的实例，可以看看，求包络的时候暂时没有考虑边界条件，自己可以完善：
fs=30;
t=0:1/fs:200;
x6=sin(2*pi*2*t)+sin(2*pi*4*t);
plot(t,x6);
xlim([0 5])
hold on
d = diff(x6);
n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;
indmax = find(d1.*d2<0 & d1>0)+1;
envmin = spline(t(indmin),x6(indmin),t);
envmax = spline(t(indmax),x6(indmax),t);
plot(t,envmin,'r');
plot(t,envmax,'m'); 参考技术B 这个函数是二维函数还是三维函数啊，若是二维直接用plot就可以画函数曲线了三维函数也有画三维的函数参考技术C 可以先通过求导找出极点,然后对极点进行函数拟合.
思路是这样.具体的代码,你懂的.....追问

呵呵具体点行不？

追答

呵呵,具体的函数我也忘记了....你好好的找一下吧

数字信号去噪基于matlab鲸鱼算法优化VMD信号去噪(目标函数为包络熵局部极小值)含Matlab源码 2091期

一、鲸鱼算法优化VMD信号去噪简介

1 VMD算法
VMD(variational mode decomposition)是一种将信号非递归分解成数个具有准正交性的带限本征模态函数的算法，其基本原理可以表示为一个约束变分问题的求解，构造的约束变分问题表示为

式中：uK为VMD分解后的第K个IMF分量；ωK为第K个IMF分量的瞬时频率；δ(t)为狄拉克函数；∗为卷积符号；

引入二次惩罚因子α和拉格朗日乘子λ使变分问题不受约束，构造的增广拉格朗日函数如下：

采用乘子交替方向算法交替更新得到K个IMF及其对应的中心频率。

从式(2)可以看出：K和α影响着VMD的分解性能。若设置的K较小，则信号的多个分量可能同时包含在1个模态中；若K较大，则会导致1个分量包含在多个模态中，迭代得到的中心频率也会重叠。对α而言，若α很大，则带宽限制就会很窄，从而导致有用的频率成分被消除；反之，冗余频率成分将会被保留下来。因此，本文提出鲸鱼优化算法来优化确定最佳参数组合(K，α)。

2 WOA算法
鲸鱼优化算法是一种基于座头鲸种群迭代进化搜索的元启发式优化算法。通过模拟座头鲸独特的气泡网攻击狩猎行为，实现复杂优化问题寻优，具体优化机理包括局部开发阶段和全局探索阶段2个部分。

1)在局部开发阶段，鲸鱼假设当前的猎物位置为目标位置，在收缩包围机制和螺旋更新机制之间以50%的概率不断靠近猎物获得最优解，其数学模型如下：

式中：X为当前的目标位置；X为鲸鱼位置；A为绝对值小于等于1的系数；D’=|X(t)-X (t)|表示鲸鱼到猎物的距离；b为常数，定义为对数螺旋的形状；l是(-1,1)中的随机数，使得螺旋形状不规则，以便更好地搜索最优解。

2)在全局探索阶段，鲸鱼以其他鲸鱼位置为参考，不断更新自身位置进行随机搜索，其数学模型可表示为

式中：A的绝对值大于1;Xrand为随机选择的鲸鱼位置。

3 基于WOA的VMD参数优化
本文提出采用鲸鱼优化算法实现变分模态分解算法的参数优化。该算法中，参数优化的关键是选择适应度函数。排列熵作为一种度量混沌时间序列复杂性的平均熵参数，时间序列复杂度参数计算更简单，抗干扰能力更强，具有更好地鲁棒性的特点，尤其适用于存在动态噪声和观测噪声的信号。基于此，本文选择排列熵作为WOA-VMD算法的适应度函数，具体构建步骤如下。

1)给定1个离散时间序列x(i),i=1～N，对该序列进行相空间重构，得到如下重构矩阵：

式中：r=N-(m-1)τ，r为重构行向量的个数；j为重构矩阵的第j行分量，j=1～r;m为嵌入维数；τ为延迟时间。

2)将重构矩阵中的r个行向量分别按照元素数值进行升序排列，得到各自的符号序列S (q)=(j1,j2,…,jm)。其中，q=1～r,r≤m!;m！为m维相空间映射的符号序列的总数；j1,j1，…，jm表示各元素在原重构分量中的索引号。

3)计算r种不同符号序列S(q)出现的概率P1,P2，…，Pr，则时间序列x(i),i=1～N的排列熵定义为

基于上述提出的适应度函数，鲸鱼优化算法和变分模态分解算法联合后的信号优化分解的流程图如图1所示，具体步骤如下。

1)输入信号，设置VMD算法中K和α的参数范围，初始化WOA模型中的各项参数，包括种群规模、最大迭代次数、空间维度以及初始种群个体。

2)对信号进行VMD分解，利用式(6)计算初始种群中每个个体的适应度。作为该方法的适应度函数，排列熵被用来衡量参数组合的分解效果。熵越小，则时间序列分布越有规律，表示VMD处理得到的IMF包含更多的有效信息；反之，时间序列越接近随机分布，IMF中噪声成分更多。因此，当排列熵最小时，对应的参数最优。

图1 WOA算法优化VMD参数的流程图
3)当|A|≤1时，选择最小排列熵对应的鲸鱼位置作为局部开发的目标值，然后根据p，选择式(3)更新鲸鱼个体的位置；当|A|>1时，随机选择一个鲸鱼的位置，根据式(4)更新鲸鱼个体的位置，保留最优适应度及对应的参数组合。

4)保留更新后的鲸鱼种群位置作为新一轮的初始种群，循环迭代，直到达到所设定的最大迭代次数为止。

5)输出最优鲸鱼个体及对应的适应度。

4 相关模态的选择
为了有效识别和分离经VMD分解后的IMF分量中的噪声分量，在对IMF分量进行时频分析的基础上引入相关系数，进一步判断IMF与原信号之间的相关程度。相关系数C的定义如下:

式中：ui(t)和x(t)为分别表示VMD分解后的第i个IMF分量和原始瞬变电磁响应电压信号；E和D表示数学上的期望和方差。

C越大，对应的IMF分量与原信号相关性越好。本文首先根据式(7)计算各个分量与原始信号的相关系数，然后将相关系数按照分量顺序进行排列，取第1个极值点位置作为有效信号分量和噪声分量的分界点。由于VMD分解得到1组从低频到高频的IMF分量，而EEMD的频率随着IMF阶数增加而逐渐减小，因此，对于VMD方法选取极值点前的分量将作为有效分量，EEMD方法极值点前对应的则为噪声分量，最后将有效分量重构得到去噪后的信号。

二、部分源代码

clc
clear
close all
warning off
format long
addpath(‘toolf’)

% 读取数据
[file,path,indx] = uigetfile(’.xlsx’;’.xls’;’.txt’;’.*’,‘File Selector’);
if indx == 1||indx==2
data=xlsread(strcat(path, file));
elseif indx == 3
data=load(strcat(path, file));
else
warning(‘请选择数据集进行实验’);
return;
end

%采样频率
fs=12800;
% 读取前1000长度的信号
len=1000;
s=data(1:len);
% 采样时间
t = (0:len-1)/fs;

%% 设定鲸鱼优化算法参数
popsize =10; % 种群大小
iter = 30; % 最大迭代次数
dim = 2; % 变量个数
lb = [100 3]; % alpha范围 K范围下限
ub = [2000 7]; % 上限

%% 鲸鱼优化算法WOA优化VMD参数

tic , % 开始计时

WOA_VMD(popsize,iter,dim,lb,ub,0); % 0表示不保存IMF，1，导出IMF并保存

toc, % 结束计时