[原创]偏微分方程数值解法的MATLAB源码【更新完毕】
说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序谢谢大家的支持！
其他的数值算法见：
..//Announce/Announce.asp?BoardID=209&id=8245004
1、古典显式格式求解抛物型偏微分方程（一维热传导方程）
function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)
%古典显式格式求解抛物型偏微分方程
%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)
%
%方程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT
%初值条件：u(x,0)=phi(x)
%边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t)
%
%输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……
% x -空间变量
% t -时间变量
%输入参数：uX -空间变量x的取值上限
% uT -时间变量t的取值上限
% phi -初值条件，定义为内联函数
% psi1 -边值条件，定义为内联函数
% psi2 -边值条件，定义为内联函数
% M -沿x轴的等分区间数
% N -沿t轴的等分区间数
% C -系数，默认情况下C=1
%
%应用举例：
%uX=1;uT=0.2;M=15;N=100;C=1;
%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');
%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);
%设置参数C的默认值
if nargin==7
C=1;
end
%计算步长
dx=uX/M;%x的步长
dt=uT/N;%t的步长
x=(0:M)*dx;
t=(0:N)*dt;
r=C*dt/dx/dx;%步长比
r1=1-2*r;
if r > 0.5
disp('r > 0.5,不稳定')
end
%计算初值和边值
U=zeros(M+1,N+1);
for i=1:M+1
U(i,1)=phi(x(i));
end
for j=1:N+1
U(1,j)=psi1(t(j));
U(M+1,j)=psi2(t(j));
end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j);
end
end
U=U';
%作出图形
mesh(x,t,U);
title('古典显式格式，一维热传导方程的解的图像') xlabel('空间变量x')
ylabel('时间变量t')
zlabel('一维热传导方程的解U')
return;
古典显式格式不稳定情况
古典显式格式稳定情况
2、古典隐式格式求解抛物型偏微分方程（一维热传导方程）
function [U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C)
%古典隐式格式求解抛物型偏微分方程
%[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C)
%
%方程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT
%初值条件：u(x,0)=phi(x)
%边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t)
%
%输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……% x -空间变量
% t -时间变量
%输入参数：uX -空间变量x的取值上限
% uT -时间变量t的取值上限
% phi -初值条件，定义为内联函数
% psi1 -边值条件，定义为内联函数
% psi2 -边值条件，定义为内联函数
% M -沿x轴的等分区间数
% N -沿t轴的等分区间数
% C -系数，默认情况下C=1
%
%应用举例：
%uX=1;uT=0.2;M=50;N=50;C=1;
%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');
%[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C);
%设置参数C的默认值
if nargin==7
C=1;
end
%计算步长
dx=uX/M;%x的步长
dt=uT/N;%t的步长
x=(0:M)*dx;
t=(0:N)*dt;
r=C*dt/dx/dx;%步长比
Diag=zeros(1,M-1);%矩阵的对角线元素
Low=zeros(1,M-2);%矩阵的下对角线元素
Up=zeros(1,M-2);%矩阵的上对角线元素
for i=1:M-2
Diag(i)=1+2*r;
Low(i)=-r;
Up(i)=-r;
end
Diag(M-1)=1+2*r;
%计算初值和边值
U=zeros(M+1,N+1);
for i=1:M+1
U(i,1)=phi(x(i));
end
for j=1:N+1
U(1,j)=psi1(t(j));
U(M+1,j)=psi2(t(j));
end
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）
for j=1:N
b1=zeros(M-1,1);
b1(1)=r*U(1,j+1);
b1(M-1)=r*U(M+1,j+1);
b=U(2:M,j)+b1;
U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);
end
U=U';
%作出图形
mesh(x,t,U);
title('古典隐式格式，一维热传导方程的解的图像')
xlabel('空间变量x')
ylabel('时间变量t')
zlabel('一维热传导方程的解U')
return;
此算法需要使用追赶法求解三对角线性方程组，这个算法在上一篇帖子中已经给出，为了方便，再给出来追赶法解三对角线性方程组
function x=EqtsForwardAndBackward(L,D,U,b)
%追赶法求解三对角线性方程组Ax=b
%x=EqtsForwardAndBackward(L,D,U,b)
%x:三对角线性方程组的解
%L:三对角矩阵的下对角线，行向量
%D:三对角矩阵的对角线，行向量
%U:三对角矩阵的上对角线，行向量
%b:线性方程组Ax=b中的b，列向量
%
%应用举例:
%L=[-1 -2 -3];D=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]';
%x=EqtsForwardAndBackward(L,D,U,b)
%检查参数的输入是否正确
n=length(D);m=length(b);
n1=length(L);n2=length(U);
if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m disp('输入参数有误！')
x=' ';
return;
end
%追的过程
for i=2:n
L(i-1)=L(i-1)/D(i-1);
D(i)=D(i)-L(i-1)*U(i-1);
end
x=zeros(n,1);
x(1)=b(1);
for i=2:n
x(i)=b(i)-L(i-1)*x(i-1);
end
%赶的过程
x(n)=x(n)/D(n);
for i=n-1:-1:1
x(i)=(x(i)-U(i)*x(i+1))/D(i);
end
return;
古典隐式格式
在以后的程序中，我们都取C=1，不再作为一个输入参数处理
3、Crank-Nicolson隐式格式求解抛物型偏微分方程
需要调用追赶法的程序
function [U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)
%Crank-Nicolson隐式格式求解抛物型偏微分方程
%[U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)
%
%方程：u_t=u_xx 0 <= x <= uX,0 <= t <= uT
%初值条件：u(x,0)=phi(x)
%边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t)
%
%输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层……% x -空间变量
% t -时间变量
%输入参数：uX -空间变量x的取值上限
% uT -时间变量t的取值上限
% phi -初值条件，定义为内联函数
% psi1 -边值条件，定义为内联函数
% psi2 -边值条件，定义为内联函数
% M -沿x轴的等分区间数
% N -沿t轴的等分区间数
%
%应用举例：
%uX=1;uT=0.2;M=50;N=50;
%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');
%[U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N);
%计算步长
dx=uX/M;%x的步长
dt=uT/N;%t的步长
x=(0:M)*dx;
t=(0:N)*dt;
r=dt/dx/dx;%步长比
Diag=zeros(1,M-1);%矩阵的对角线元素
Low=zeros(1,M-2);%矩阵的下对角线元素
Up=zeros(1,M-2);%矩阵的上对角线元素
for i=1:M-2
Diag(i)=1+r;
Low(i)=-r/2;
Up(i)=-r/2;
end
Diag(M-1)=1+r;
%计算初值和边值
U=zeros(M+1,N+1);
for i=1:M+1
U(i,1)=phi(x(i));
end
for j=1:N+1
U(1,j)=psi1(t(j));
U(M+1,j)=psi2(t(j));
end
B=zeros(M-1,M-1);
for i=1:M-2
B(i,i)=1-r;
B(i,i+1)=r/2;
B(i+1,i)=r/2;
end
B(M-1,M-1)=1-r;
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）for j=1:N
b1=zeros(M-1,1);
b1(1)=r*(U(1,j+1)+U(1,j))/2;
b1(M-1)=r*(U(M+1,j+1)+U(M+1,j))/2;
b=B*U(2:M,j)+b1;
U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);
end
U=U';
%作出图形
mesh(x,t,U);
title('Crank-Nicolson隐式格式，一维热传导方程的解的图像')
xlabel('空间变量x')
ylabel('时间变量t')
zlabel('一维热传导方程的解U')
return;
Crank-Nicolson隐式格式
4、正方形区域Laplace方程Diriclet问题的求解
需要调用Jacobi迭代法和Guass-Seidel迭代法求解线性方程组
function [U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type) %正方形区域Laplace方程的Diriclet边值问题的差分求解
%此程序需要调用Jacobi迭代法或者Guass-Seidel迭代法求解线性方程组
%[U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)
%
%方程：u_xx+u_yy=0 0<=x,y<=ub
%边值条件：u(0,y)=phi1(y)
% u(ub,y)=phi2(y)
% u(x,0)=psi1(x)
% u(x,ub)=psi2(x)
%
%输出参数：U -解矩阵，第一行表示y=0时的值，第二行表示第y=h时的值……% x -横坐标
% y -纵坐标
%输入参数：ub -变量边界值的上限
% phi1,phi2,psi1,psi2 -边界函数，定义为内联函数
% M -横纵坐标的等分区间数
% type -求解差分方程的迭代格式，若type='Jacobi'，采用Jacobi迭代格式% 若type='GS'，采用Guass-Seidel迭代格式。默认情况下，type='GS' %
%应用举例：
%ub=4;M=20;
%phi1=inline('y*(4-y)');phi2=inline('0');psi1=inline('sin(pi*x/4)');psi2=inline('0');
%[U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,'GS');
if nargin==6
type='GS';
end
%步长
h=ub/M;
%横纵坐标
x=(0:M)*h;
y=(0:M)*h;
%差分格式的矩阵形式AU=K
%构造矩阵A
M2=(M-1)^2;
A=zeros(M2);
for i=1:M2
A(i,i)=4;
end
for i=1:M2-1
if mod(i,M-1)~=0
A(i,i+1)=-1;
A(i+1,i)=-1;
end
end
for i=1:M2-M+1
A(i,i+M-1)=-1;
A(i+M-1,i)=-1;
end
U=zeros(M+1);
%边值条件
for i=1:M+1
U(i,1)=psi1((i-1)*h);
U(i,M+1)=psi2((i-1)*h);
U(1,i)=phi1((i-1)*h);
U(M+1,i)=phi2((i-1)*h);
end
%构造K
K=zeros(M2,1);
for i=1:M-1
K(i)=U(i+1,1);
K(M2-i+1)=U(i+1,M+1);
end
K(1)=K(1)+U(1,2);
K(M-1)=K(M-1)+U(M+1,2);
K(M2-M+2)=K(M2-M+2)+U(1,M);
K(M2)=K(M2)+U(M+1,M);
for i=2:M-2
K((M-1)*(i-1)+1)=U(1,i+1);
K((M-1)*i)=U(M+1,i+1);
end
x0=ones(M2,1);
switch type
%调用Guass-Seidel迭代法求解线性方程组AU=K case 'Jacobi'
X=EqtsJacobi(A,K,x0);
%调用Guass-Seidel迭代法求解线性方程组AU=K case 'GS'
X=EqtsGS(A,K,x0);
otherwise
disp('差分格式类型输入错误')
return;
end
%把求解结果化成矩阵型式
for i=2:M
for j=2:M
U(j,i)=X(j-1+(M-1)*(i-2));
end
end
U=U';
%作出图形
mesh(x,y,U);
title('五点差分格式Laplace方程Diriclet问题的解的图像')
xlabel('x')
ylabel('y')
zlabel('Laplace方程Diriclet问题的解U')
return;
正方形区域Laplace方程五点差分格式
5、一阶双曲型方程的差分方法
function [U x t]=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)
%一阶双曲型方程的差分格式
%[U x t]=PDEHyperbolic(uX,uT,M,N,C,phi,psi1,psi2,type)
%
%方程:u_t+C*u_x=0 0 <= t <= uT, 0 <= x <= uX
%初值条件:u(x,0)=phi(x)
%
%输出参数:U -解矩阵，第一行表示初值，第二行表示第2个时间层……
% x -横坐标
% t -纵坐标，时间
%输入参数:uX -变量x的上界
% uT -变量t的上界
% M -变量x的等分区间数
% N -变量t的等分区间数
% C -系数
% phi -初值条件函数，定义为内联函数
% psi1,psi2 -边值条件函数，定义为内联函数
% type -差分格式，从下列值中选取
% -type='LaxFriedrichs'，采用Lax-Friedrichs差分格式求解
% -type='CourantIsaacsonRees'，采用Courant-Isaacson-Rees差分格式求解
% -type='LeapFrog'，采用Leap-Frog（蛙跳）差分格式求解
% -type='LaxWendroff'，采用Lax-Wendroff差分格式求解
% -type='CrankNicolson'，采用Crank-Nicolson差分格式求解，此格式需调用追赶法% 求解三对角线性方程组
%
h=uX/M;%变量x的步长
k=uT/N;%变量t的步长
r=k/h;%步长比
x=(0:M)*h;
t=(0:N)*k;
U=zeros(M+1,N+1);
%初值条件
for i=1:M+1
U(i,1)=phi(x(i));
end
%边值条件
for j=1:N+1
U(1,j)=psi1(t(j));
U(M+1,j)=psi2(t(j));
%U(1,j)=NaN;
%U(M+1,j)=NaN;
end
switch type
%Lax-Friedrichs差分格式
case 'LaxFriedrichs'
if abs(C*r)>1
disp('|C*r|>1，Lax-Friedrichs差分格式不稳定！')
end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=(U(i+1,j)+U(i-1,j))/2-C*r*(U(i+1,j)-U(i-1,j))/2;
end
end
%Courant-Isaacson-Rees差分格式
case 'CourantIsaacsonRees'
if C<0
disp('C<0，采用前差公式')
if C*r<-1
disp('Courant-Isaacson-Lees差分格式不稳定！') end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=(1+C*r)*U(i,j)-C*r*U(i+1,j);
end
end
else
disp('C>0，采用后差公式')
if C*r>1
disp('Courant-Isaacson-Lees差分格式不稳定！') end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=C*r*U(i-1,j)+(1-C*r)*U(i,j);
end
end
end
%Leap-Frog（蛙跳）差分格式
case 'LeapFrog'
phi2=input('请输入第二层初值条件函数：psi2=');
if abs(C*r)>1
disp('|C*r|>1，Leap-Frog差分格式不稳定！')
end
%第二层初值条件
for i=1:M+1
U(i,2)=phi2(x(i));
end
%逐层求解
for j=2:N
for i=2:M
U(i,j+1)=U(i,j-1)-C*r*(U(i+1,j)-U(i-1,j));
end
end
%Lax-Wendroff差分格式
case 'LaxWendroff'
if abs(C*r)>1
disp('|C*r|>1，Lax-Wendroff差分格式不稳定！')
end
%逐层求解
for j=1:N
for i=2:M
U(i,j+1)=U(i,j)-C*r*(U(i+1,j)-U(i-1,j))/2+C^2*r^2*(U(i+1,j)-2*U(i,j)+U(i-1,j))/2;
end
end
%Crank-Nicolson隐式差分格式，需调用追赶法求解三对角线性方程组的算法
case 'CrankNicolson'
Diag=zeros(1,M-1);%矩阵的对角线元素
Low=zeros(1,M-2);%矩阵的下对角线元素
Up=zeros(1,M-2);%矩阵的上对角线元素
for i=1:M-2
Diag(i)=4;
Low(i)=-r*C;
Up(i)=r*C;
end
Diag(M-1)=4;
B=zeros(M-1,M-1);
for i=1:M-2
B(i,i)=4;
B(i,i+1)=-r*C;
B(i+1,i)=r*C;
end
B(M-1,M-1)=4;
%逐层求解，需要使用追赶法（调用函数EqtsForwardAndBackward）
for j=1:N
b1=zeros(M-1,1);
b1(1)=r*C*(U(1,j+1)+U(1,j))/2;
b1(M-1)=-r*C*(U(M+1,j+1)+U(M+1,j))/2;
b=B*U(2:M,j)+b1;
U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);
end
otherwise
disp('差分格式类型输入有误！')
return;
end
U=U';
%作出图形
mesh(x,t,U);
title([type '格式求解一阶双曲型方程的解的图像']);
xlabel('空间变量x');
ylabel('时间变量t');
zlabel('一阶双曲型方程的解U');
return;
1.硬币模拟试验 源代码： clear; clc; head_count=0; p1_hist= [0]; p2_hist= [0]; n = 1000; p1 = 0.3; p2=0.03; head = figure(1); rand('seed',sum(100*clock)); fori = 1:n tmp = rand(1); if(tmp<= p1) head_count = head_count + 1; end p1_hist (i) = head_count /i; end figure(head); subplot(2,1,1); plot(p1_hist); grid on; hold on; xlabel('重复试验次数'); ylabel('正面向上的比率'); title('p=0.3试验次数N与正面向上比率的函数图'); head_count=0; fori = 1:n tmp = rand(1); if(tmp<= p2) head_count = head_count + 1; end p2_hist (i) = head_count /i; end figure(head); subplot(2,1,2); plot(p2_hist); grid on; hold on; xlabel('重复试验次数'); ylabel('正面向上的比率'); title('p=0.03试验次数N与正面向上比率的函数图'); 实验结果：
2.不同次数的随机试验均值方差比较 源代码： clear ; clc; close; rand('seed',sum(100*clock)); Titles = ['n=5时' 'n=20时' 'n=25时' 'n=50时' 'n=100时']; Titlestr = cellstr(Titles); X_n_bar=[0]; %the samples of the X_n_bar X_n=[0]; %the samples of X_n N=[5,10,25,50,100]; j=1; num_X_n = 100; num_X_n_bar = 100; h_X_n_bar = figure(1);
数字图像处理MATLAB函数源代码MATLAB函数源代码 function varargout = DIP(varargin) % DIP MATLAB code for DIP.fig % DIP, by itself, creates a new DIP or raises the existing % singleton*. % % Edit the above text to modify the response to help DIP % Last Modified by GUIDE v2.5 27-May-2012 11:43:05 gui_Singleton = 1; gui_State = struct('gui_Name', mfilename, ... 'gui_Singleton', gui_Singleton, ... 'gui_OpeningFcn', @DIP_OpeningFcn, ... 'gui_OutputFcn', @DIP_OutputFcn, ... 'gui_LayoutFcn', [] , ... 'gui_Callback', []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end
例错误！文档中没有指定样式的文字。-1 %周期信号（方波）的展开,fb_jinshi.m close all; clear all; N=100; %取展开式的项数为2N＋1项 T=1; fs=1/T; N_sample=128; %为了画出波形，设置每个周期的采样点数 dt = T/N_sample; t=0:dt:10*T-dt; n=-N:N; Fn = sinc(n/2).*exp(-j*n*pi/2); Fn(N+1)=0; ft = zeros(1,length(t)); for m=-N:N ft = ft + Fn(m+N+1)*exp(j*2*pi*m*fs*t); end plot(t,ft) 例错误！文档中没有指定样式的文字。-4 利用FFT计算信号的频谱并与信号的真实频谱的抽样比较。 脚本文件T2F.m定义了函数T2F，计算信号的傅立叶变换。 function [f,sf]= T2F(t,st) %This is a function using the FFT function to calculate a signal's Fourier %Translation %Input is the time and the signal vectors,the length of time must greater %than 2 %Output is the frequency and the signal spectrum dt = t(2)-t(1); T=t(end); df = 1/T; N = length(st); f=-N/2*df:df:N/2*df-df; sf = fft(st); sf = T/N*fftshift(sf); 脚本文件F2T.m定义了函数F2T，计算信号的反傅立叶变换。 function [t st]=F2T(f,sf) %This function calculate the time signal using ifft function for the input %signal's spectrum
1.绘制云图 Ex=18 En=2 He=0.2 hold on for i=1:1000 Enn=randn(1)*He+En; x(i)=randn(1)*Enn+Ex; y(i)=exp(-(x(i)-Ex)^2/(2*Enn^2)); plot(x(i),y(i),'*') end Ex=48.7 En=9.1 He=0.39 hold on for i=1:1000 Enn=randn(1)*He+En; x(i)=randn(1)*Enn+Ex; y(i)=exp(-(x(i)-Ex)^2/(2*Enn^2)); plot(x(i),y(i),'*')
end 2.求期望、熵及超熵 X1=[51.93 52.51 54.70 43.14 43.85 44.48 44.61 52.08]; Y1=[0.91169241573 0.921875 0.96032303371 0.75737359551 0.76983848315 0.7808988764 0.78318117978 0.9143258427]; m=8; Ex=mean(X1) En1=zeros(1,m); for i=1:m En1(1,i)=abs(X1(1,i)-Ex)/sqrt(-2*log(Y1(1,i))); end En=mean(En1); He=0; for i=1:m He=He+(En1(1,i)-En)^2; end En=mean(En1) He=sqrt(He/(m-1)) 3.平顶山so2环境: X1=[0.013 0.04 0.054 0.065 0.07 0.067 0.058 0.055 0.045]; Y1=[0.175675676 0.540540541 0.72972973 0.878378378
matlab语音识别系统(源代码)最新版
目录 一、设计任务及要求 (1) 二、语音识别的简单介绍 2.1语者识别的概念 (2) 2.2特征参数的提取 (3) 2.3用矢量量化聚类法生成码本 (3) 2.4VQ的说话人识别 (4) 三、算法程序分析 3.1函数关系 (4) 3.2代码说明 (5) 3.2.1函数mfcc (5) 3.2.2函数disteu (5) 3.2.3函数vqlbg (6) 3.2.4函数test (6) 3.2.5函数testDB (7) 3.2.6 函数train (8) 3.2.7函数melfb (8) 四、演示分析 (9) 五、心得体会 (11) 附：GUI程序代码 (12)
一、设计任务及要求 用MATLAB实现简单的语音识别功能； 具体设计要求如下： 用MATLAB实现简单的数字1~9的语音识别功能。 二、语音识别的简单介绍 基于VQ的说话人识别系统，矢量量化起着双重作用。在训练阶段，把每一个说话者所提取的特征参数进行分类，产生不同码字所组成的码本。在识别(匹配)阶段，我们用VQ方法计算平均失真测度(本系统在计算距离d时，采用欧氏距离测度)，从而判断说话人是谁。 语音识别系统结构框图如图1所示。 图1 语音识别系统结构框图 2.1语者识别的概念 语者识别就是根据说话人的语音信号来判别说话人的身份。语音是人的自然属性之一，由于说话人发音器官的生理差异以及后天形成的行为差异，每个人的语音都带有强烈的个人色彩，这就使得通过分析语音信号来识别说话人成为可能。用语音来鉴别说话人的身份有着许多独特的优点，如语音是人的固有的特征，不会丢失或遗忘；语音信号的采集方便，系统设备成本低；利用电话网络还可实现远程客户服务等。因此，近几年来，说话人识别越来越多的受到人们的重视。与其他生物识别技术如指纹识别、手形识别等相比较，说话人识别不仅使用方便，而且属于非接触性，容易被用户接受，并且在已有的各种生物特征识别技术中，是唯一可以用作远程验证的识别技术。因此，说话人识别的应用前景非常广泛：今天，说话人识别技术已经关系到多学科的研究领域，不同领域中的进步都对说话人识别的发展做出了贡献。说话人识别技术是集声学、语言学、计算机、信息处理和人工智能等诸多领域的一项综合技术，应用需求将十分广阔。在吃力语音信号的时候如何提取信号中关键的成分尤为重要。语音信号的特征参数的好坏直接导致了辨别的准确性。
% *-*--*-*-*-*-*-*-*-*-*-*-*图像处理*-*-*-*-*-*-*-*-*-*-*-* %{ % （一）图像文件的读/写 A=imread('drum.jpg'); % 读入图像 imshow(A); % 显示图像 imwrite(A,'drum.jpg'); info=imfinfo('drum.jpg') % 查询图像文件信息 % 用colorbar函数将颜色条添加到坐标轴对象中 RGB=imread('drum.jpg'); I=rgb2gray(RGB); % 把RGB图像转换成灰度图像 h=[1 2 1;0 0 0;-1 -2 -1]; I2=filter2(h,I); imshow(I2,[]); colorbar('vert') % 将颜色条添加到坐标轴对象中 % wrap函数将图像作为纹理进行映射 A=imread('4.jpg'); imshow(A); I=rgb2gray(RGB); [x,y,z]=sphere; warp(x,y,z,I); % 用warp函数将图像作为纹理进行映射 %} % subimage函数实现一个图形窗口中显示多幅图像 RGB=imread('drum.jpg'); I=rgb2gray(RGB); subplot(1,2,1); subimage(RGB); % subimage函数实现一个图形窗口中显示多幅图像subplot(1,2,2),subimage(I); % *-*--*-*-*-*-*-*-*-*-*-*-*图像处理*-*-*-*-*-*-*-*-*-*-*-* % （二）图像处理的基本操作 % ----------------图像代数运算------------------ %{ % imadd函数实现两幅图像的相加或给一幅图像加上一个常数 % 给图像每个像素都增加亮度 I=imread('4.jpg'); J=imadd(I,100); % 给图像增加亮度 subplot(1,2,1),imshow(I);title('原图'); subplot(1,2,2),imshow(J);title('增加亮度图'); % % imsubtract函数实现将一幅图像从另一个图像中减去或减去一个常数I=imread('drum.jpg'); J=imsubtract(I,100); % 给图像减去亮度 subplot(1,2,1),imshow(I);
模型((P C2R)的MATLAB程序 clear X=[]; %用户输入多指标输入矩阵X Y=[]; %用户输入多指标输出矩阵Y n=size(X',1); m=size(X,1); s=size(Y,1); A=[-X' Y']; b=zeros(n, 1); LB=zeros(m+s,1); UB=[]; for i=1:n; f= [zeros(1,m) -Y(:,i)']; Aeq=[X(:,i)' zeros(1,s)]; beq=1; w(:,i)=LINPROG(f,A,b,Aeq,beq,LB,UB); %解线性规划，得DMU;的最佳权向量w; E(i, i)=Y(:,i)'*w(m+1:m+s,i); %求出DMU i的相对效率值E ii end w %输出最佳权向量 E %输出相对效率值E ii Omega=w(1:m,:) %输出投入权向量。 mu=w(m+1:m+s,:) %输出产出权向量。 模型(D C2R)的MATLAB程序 clear X=[]; %用户输入多指标输入矩阵X Y=[]; %用户输入多指标输出矩阵Y n=size(X',1); m=size(X,1); s=size(Y,1); epsilon=10^-10; %定义非阿基米德无穷小 =10-10 f=[zeros(1,n) -epsilon*ones(1,m+s) 1]; %目标函数的系数矩阵： 的系数为0,s-,s+的系数为- e, 的系数为1； A=zeros(1,n+m+s+1); b=0; %<=约束； LB=zeros(n+m+s+1,1); UB=[]; %变量约束； LB(n+m+s+1)= -Inf; %-Inf表示下限为负无穷大。 for i=1:n; Aeq=[X eye(m) zeros(m,s) -X(:,i) Y zeros(s,m) -eye(s) zeros(s,1)]; beq=[zeros(m, 1 ) Y(:,i)]; w(:,i)=LINPROG (f,A,b,Aeq,beq,LB,UB); %解线性规划，得DMU的最佳权向量w; end w %输出最佳权向量 lambda=w(1:n,:) %输出 s_minus=w(n+1:n+m,:) %输出s- s_plus=w(n+m+1:n+m+s,:) %输出s+ theta=w(n+m+s+1,:) %输出
%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算********** clear clc load E:\data\IEEE014_Node.txt Node=IEEE014_Node; weishu=size(Node); nnum=weishu(1,1); %节点总数 load E:\data\IEEE014_Branch.txt branch=IEEE014_Branch; bwei=size(branch); bnum=bwei(1,1); %支路总数 Y=(zeros(nnum)); Sj=100; %********************************节点导纳矩阵******************************* for m=1:bnum; s=branch(m,1); %首节点 e=branch(m,2); %末节点 R=branch(m,3); %支路电阻 X=branch(m,4); %支路电抗 B=branch(m,5); %支路对地电纳 k=branch(m,6); if k==0 %无变压器支路情形 Y(s,e)=-1/(R+j*X); %互导纳 Y(e,s)=Y(s,e); end if k~=0 %有变压器支路情形 Y(s,e)=-(1/((R+j*X)*k)); Y(e,s)=Y(s,e); Y(s,s)=-(1-k)/((R+j*X)*k^2); Y(e,e)=-(k-1)/((R+j*X)*k); %对地导纳 end Y(s,s)=Y(s,s)-j*B/2; Y(e,e)=Y(e,e)-j*B/2; %自导纳的计算情形 end for t=1:nnum; Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13); %求支路自导纳 end G=real(Y); %电导 B=imag(Y); %电纳 %******************节点分类************************************* * pq=0; pv=0; blancenode=0; pqnode=zeros(1,nnum); pvnode=zeros(1,nnum); for m=1:nnum; if Node(m,2)==3 blancenode=m; %平衡节点编号 else if Node(m,2)==0 pq=pq+1; pqnode(1,pq)=m; %PQ 节点编号 else if Node(m,2)==2 pv=pv+1; pvnode(1,pv)=m; %PV 节点编号 end end end end %*****************************设置电压初值********************************** Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化 for n=1:nnum Uoriginal(1,n)=Node(n,9); %对各点电压赋初值 if Node(n,9)==0;
数字图像处理实验一 15生医 一、实验内容 产生右图所示图像 f1(m,n)，其中图像大小为256 ×256，中间亮条为128×32，暗处=0，亮处=100。 对其进行FFT： ①同屏显示原图f1(m,n)和FFT(f1)的幅度谱图； ②若令f2(m,n)=（-1）^（m+n）f1(m,n)，重复 以上过程，比较二者幅度谱的异同，简述理由； ③若将f2(m,n)顺时针旋转90度得到f3(m,n)，试显示FFT(f3)的 幅度谱，并与FFT(f2)的幅度谱进行比较； ④若将f1(m,n) 顺时针旋转90度得到f4(m,n)，令f5(m,n) = f1(m,n) + f4(m,n)，试显示FFT(f5)的幅度谱，指出其与 FFT(f1)和FFT(f4)的关系； ⑤若令f6(m,n)=f2(m,n)+f3(m,n)，试显示FFT(f6)的幅度谱，并指出其与 FFT(f2)和FFT(f3)的关系，比较FFT(f6)和FFT(f5)的幅度谱。 二、运行环境 MATLAB R2014a 三、运行结果及分析 1.同屏显示原图f1(m,n)和FFT(f1)的幅度谱图：
50100150200250 100150200250 50100150200250 100150200250 2.令f2(m,n)=（-1）^（m+n ）f1(m,n)，对其进行FFT ，比较f2与f1幅度谱的异同，简述理由： 50100150200250 100150200250 50100150200250 100150200250 异同及理由：①空域：f2由于前边乘了系数（-1）^（m+n ），导致灰度值有正有负，而在MATLAB 的imshow 函数中默认把负值变为0（有些情况是取反），所以形成了如左图所示的黑白花纹。②频域：FFT(2)
附录 ＭＡＴＬＡＢ图像处理命令 １．ａｐｐｌｙｌｕｔ 功能： 在二进制图像中利用lookup表进行边沿操作。 语法： A = applylut(BW,lut) 举例 lut = makelut('sum(x(:)) == 4',2); BW1 = imread('text.tif'); BW2 = applylut(BW1,lut); imshow(BW1) figure, imshow(BW2) 相关命令： makelut ２．ｂｅｓｔｂｌｋ 功能： 确定进行块操作的块大小。 语法： siz = bestblk([m n],k) [mb,nb] = bestblk([m n],k) 举例 siz = bestblk([640 800],72) siz = 64 50 相关命令： blkproc ３．ｂｌｋｐｒｏｃ 功能：
MATLAB 高级应用——图形及影像处理 320 实现图像的显式块操作。 语法： B = blkproc(A,[m n],fun) B = blkproc(A,[m n],fun,P1,P2,...) B = blkproc(A,[m n],[mborder nborder],fun,...) B = blkproc(A,'indexed',...) 举例 I = imread('alumgrns.tif'); I2 = blkproc(I,[8 8],'std2(x)*ones(size(x))'); imshow(I) figure, imshow(I2,[]); 相关命令： colfilt, nlfilter,inline ４．ｂｒｉｇｈｔｅｎ 功能： 增加或降低颜色映像表的亮度。 语法： brighten(beta) newmap = brighten(beta) newmap = brighten(map,beta) brighten(fig,beta) 相关命令： imadjust, rgbplot ５．ｂｗａｒｅａ 功能： 计算二进制图像对象的面积。 语法： total = bwarea(BW) 举例 BW = imread('circles.tif'); imshow(BW);
云模型简介及个人理解m a t l a b程序 集团档案编码：[YTTR-YTPT28-YTNTL98-UYTYNN08]
随着不确定性研究的深入，越来越多的科学家相信，不确定性是这个世界的魅力所在，只有不确定性本身才是确定的。在众多的不确定性中，和是最基本的。针对和在处理不确定性方面的不足，1995年我国工程院院士教授在概率论和模糊数学的基础上提出了云的概念，并研究了模糊性和随机性及两者之间的关联性。自李德毅院士等人提出云模型至今，云模型已成功的应用到、、、智能控制、等众多领域. 设是一个普通集合。 , 称为论域。关于论域中的模糊集合，是指对于任意元素都存在一个有稳定倾向的随机数，叫做对的隶属度。如果论域中的元素是简单有序的，则可以看作是基础变量，隶属度在上的分布叫做隶属云；如果论域中的元素不是简单有序的，而根据某个法则，可将映射到另一个有序的论域上，中的一个且只有一个和对应，则为基础变量，隶属度在上的分布叫做隶属云[1] 。 数字特征
云模型表示自然语言中的基元——语言值，用云的数字特征——期望Ex，熵En和超熵He表示语言值的数学性质 [3] 。 期望 Ex：云滴在论域空间分布的期望，是最能够代表定性概念的点，是这个概念量化的最典型样本。 熵 En：“熵”这一概念最初是作为描述热力学的一个状态参量，此后又被引入统计物理学、信息论、复杂系统等，用以度量不确定的程度。在云模型中，熵代表定性概念的可度量粒度，熵越大，通常概念越宏观，也是定性概念不确定性的度量，由概念的随机性和模糊性共同决定。一方面, En是定性概念随机性的度量，反映了能够代表这个定性概念的云滴的离散程度；另一方面，又是定性概念亦此亦彼性的度量，反映了在论域空间可被概念接受的云滴的取值范围。用同一个数字特征来反映随机性和模糊性，也必然反映他们之间的关联性。 超熵 He：熵的不确定性度量，即熵的熵，由熵的随机性和模糊性共同决定。反映了每个数值隶属这个语言值程度的凝聚性，即云滴的凝聚程度。超熵越大，云的离散程度越大，隶属度的随机性也随之增大，云的厚度也越大。
最常用的一些图像处理Matlab源代 码 #1：数字图像矩阵数据的显示及其傅立叶变换 #2：二维离散余弦变换的图像压缩 #3：采用灰度变换的方法增强图像的对比度 #4：直方图均匀化 #5：模拟图像受高斯白噪声和椒盐噪声的影响 #6：采用二维中值滤波函数medfilt2对受椒盐噪声干扰的图像滤波 #7：采用MATLAB中的函数filter2对受噪声干扰的图像进行均值滤波 #8：图像的自适应魏纳滤波 #9：运用5种不同的梯度增强法进行图像锐化 #10：图像的高通滤波和掩模处理 #11：利用巴特沃斯（Butterworth）低通滤波器对受噪声干扰的图像进行平滑处理 #12：利用巴特沃斯（Butterworth）高通滤波器对受噪声干扰的图像进行平滑处理 1.数字图像矩阵数据的显示及其傅立叶变换 f=zeros(30,30); f(5:24,13:17)=1; imshow(f, 'notruesize'); F=fft2(f,256,256); % 快速傅立叶变换算法只能处矩阵维数为2的幂次，f矩阵不 % 是，通过对f矩阵进行零填充来调整 F2=fftshift(F); % 一般在计算图形函数的傅立叶变换时，坐标原点在 % 函数图形的中心位置处，而计算机在对图像执行傅立叶变换 % 时是以图像的左上角为坐标原点。所以使用函数fftshift进 %行修正，使变换后的直流分量位于图形的中心； figure,imshow(log(abs(F2)),[-1 5],'notruesize');
2 二维离散余弦变换的图像压缩I=imread('cameraman.tif'); % MATLAB自带的图像imshow(I); clear;close all I=imread('cameraman.tif'); imshow(I); I=im2double(I); T=dctmtx(8); B=blkproc(I,[8 8], 'P1*x*P2',T,T'); Mask=[1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; B2=blkproc(B,[8 8],'P1.*x',Mask); % 此处为点乘(.*) I2=blkproc(B2,[8 8], 'P1*x*P2',T',T); figure,imshow(I2); % 重建后的图像 3.采用灰度变换的方法增强图像的对比度I=imread('rice.tif'); imshow(I); figure,imhist(I); J=imadjust(I,[0.15 0.9], [0 1]); figure,imshow(J); figure,imhist(J);
一、编写程序完成不同滤波器的图像频域降噪和边缘增强的算法并进行比较，得出结论。 1、不同滤波器的频域降噪 1.1 理想低通滤波器(ILPF) I1=imread('eight.tif'); %读取图像 I2=im2double(I1); I3=imnoise(I2,'gaussian',0.01); I4=imnoise(I3,'salt & pepper',0.01); figure,subplot(1,3,1); imshow(I2) %显示灰度图像 title('原始图像'); %为图像添加标题 subplot(1,3,2); imshow(I4) %加入混合躁声后显示图像 title('加噪后的图像'); s=fftshift(fft2(I4)); %将灰度图像的二维不连续Fourier 变换的零频率成分 移到频谱的中心 [M,N]=size(s); %分别返回s的行数到M中，列数到N中n1=floor(M/2); %对M/2进行取整 n2=floor(N/2); %对N/2进行取整 d0=40; %初始化d0 for i=1:M for j=1:N d=sqrt((i-n1)^2+(j-n2)^2); %点（i,j）到傅立叶变换中心的距离 if d<=d0 %点（i,j）在通带内的情况 h=1; %通带变换函数 else %点（i,j）在阻带内的情况 h=0; %阻带变换函数 end s(i,j)=h*s(i,j); %ILPF滤波后的频域表示
end end s=ifftshift(s); %对s进行反FFT移动 s=im2uint8(real(ifft2(s))); %对s进行二维反离散的Fourier变换后，取复 数的实部转化为无符号8位整数 subplot(1,3,3); %创建图形图像对象 imshow(s); %显示ILPF滤波后的图像 title('ILPF滤波后的图像（d=40）'); 运行结果： 1.2 二阶巴特沃斯低通滤波器(BLPF) I1=imread('eight.tif'); %读取图像 I2=im2double(I1); I3=imnoise(I2,'gaussian',0.01); I4=imnoise(I3,'salt & pepper',0.01); figure,subplot(1,3,1); imshow(I2) %显示灰度图像 title('原始图像'); %为图像添加标题 subplot(1,3,2); imshow(I4) %加入混合躁声后显示图像 title('加噪后的图像'); s=fftshift(fft2(I4));%将灰度图像的二维不连续Fourier 变换的零频率成分 移到频谱的中心 [M,N]=size(s); %分别返回s的行数到M中，列数到N中n=2; %对n赋初值
主函数源程序（main.m） %------基本粒子群优化算法（Particle Swarm Optimization）----------- %------名称：基本粒子群优化算法（PSO） %------作用：求解优化问题 %------说明：全局性，并行性，高效的群体智能算法 %------初始格式化-------------------------------------------------- clear all; clc; format long; %------给定初始化条件---------------------------------------------- c1=1.4962; %学习因子1 c2=1.4962; %学习因子2 w=0.7298; %惯性权重 MaxDT=1000; %最大迭代次数 D=10; %搜索空间维数（未知数个数） N=40; %初始化群体个体数目 eps=10^(-6); %设置精度(在已知最小值时候用) %------初始化种群的个体(可以在这里限定位置和速度的范围)------------ for i=1:N for j=1:D x(i,j)=randn; %随机初始化位置 v(i,j)=randn; %随机初始化速度 end end %------先计算各个粒子的适应度，并初始化Pi和Pg---------------------- for i=1:N p(i)=fitness(x(i,:),D); y(i,:)=x(i,:); end pg=x(1,:); %Pg为全局最优 for i=2:N if fitness(x(i,:),D)
close all clear echo on clc % NEWFF——生成一个新的前向神经网络 % TRAIN——对 BP 神经网络进行训练 % SIM——对 BP 神经网络进行仿真 % 定义训练样本 % P为输入矢量 P=[0.7317 0.6790 0.5710 0.5673 0.5948;0.6790 0.5710 0.5673 0.5948 0.6292; ... 0.5710 0.5673 0.5948 0.6292 0.6488;0.5673 0.5948 0.6292 0.6488 0.6130; ... 0.5948 0.6292 0.6488 0.6130 0.5654; 0.6292 0.6488 0.6130 0.5654 0.5567; ... 0.6488 0.6130 0.5654 0.5567 0.5673;0.6130 0.5654 0.5567 0.5673 0.5976; ... 0.5654 0.5567 0.5673 0.5976 0.6269;0.5567 0.5673 0.5976 0.6269 0.6274; ... 0.5673 0.5976 0.6269 0.6274 0.6301;0.5976 0.6269 0.6274 0.6301 0.5803; ... 0.6269 0.6274 0.6301 0.5803 0.6668;0.6274 0.6301 0.5803 0.6668 0.6896; ... 0.6301 0.5803 0.6668 0.6896 0.7497]; % T为目标矢量 T=[0.6292 0.6488 0.6130 0.5654 0.5567 0.5673 0.5976 ... 0.6269 0.6274 0.6301 0.5803 0.6668 0.6896 0.7497 0.8094]; % Ptest为测试输入矢量 Ptest=[0.5803 0.6668 0.6896 0.7497 0.8094;0.6668 0.6896 0.7497 0.8094 0.8722; ... 0.6896 0.7497 0.8094 0.8722 0.9096]; % Ttest为测试目标矢量 Ttest=[0.8722 0.9096 1.0000]; % 创建一个新的前向神经网络 net=newff(minmax(P'),[12,1],{'logsig','purelin'},'traingdm'); % 设置训练参数 net.trainParam.show = 50; net.trainParam.lr = 0.05; net.trainParam.mc = 0.9; net.trainParam.epochs = 5000; net.trainParam.goal = 0.001; % 调用TRAINGDM算法训练 BP 网络 [net,tr]=train(net,P',T); % 对BP网络进行仿真 A=sim(net,P'); figure; plot((1993:2007),T,'-*',(1993:2007),A,'-o'); title('网络的实际输出和仿真输出结果，*为真实值，o为预测值'); xlabel('年份'); ylabel('客运量'); % 对BP网络进行测试 A1=sim(net,Ptest');
信号与系统实验报告——图像处理 学院：信息科学与工程学院 专业：2014级通信工程 组长：** 组员：** 2017.01.02
目录 目录 (2) 实验一图像一的细胞计数 (3) 一、实验内容及步骤 (3) 二、Matlab程序代码 (3) 三、数据及结果 (4) 实验二图像二的图形结构提取 (5) 一、实验内容及步骤 (5) 二、Matlab程序代码 (5) 三、数据及结果 (6) 实验三图像三的图形结构提取 (7) 一、实验内容及步骤 (7) 二、Matlab程序代码 (7) 三、数据及结果 (8) 实验四图像四的傅里叶变化及巴特沃斯低通滤波 (9) 一、实验内容及步骤 (9) 二、Matlab程序代码 (9) 三、数据及结果 (10) 实验五图像五的空间域滤波与频域滤波 (11) 一、实验内容及步骤 (11) 二、Matlab程序代码 (11) 三、数据及结果 (12)
实验一图像一的细胞计数 一、实验内容及步骤 将该图形进行一系列处理，计算得到途中清晰可见细胞的个数。 首先，由于原图为RGB三色图像处理起来较为麻烦，所以转为灰度图，再进行二值化化为黑白图像，得到二值化图像之后进行中值滤波得到细胞分布的初步图像，为了方便计数对图像取反，这时进行一次计数，发现得到的个数远远多于实际个数，这时在进行一次中值滤波，去掉一些不清晰的像素点，剩下的应该为较为清晰的细胞个数，再次计数得到大致结果。 二、Matlab程序代码 clear;close all; Image = imread('1.jpg'); figure,imshow(Image),title('原图'); Image=rgb2gray(Image); figure,imshow(Image),title('灰度图'); Theshold = graythresh(Image); Image_BW = im2bw(Image,Theshold); Reverse_Image_BW22=~Image_BW; figure,imshow(Image_BW),title('二值化图像'); Image_BW_medfilt= medfilt2(Image_BW,[3 3]); figure,imshow(Image_BW_medfilt),title('中值滤波后的二值化图像'); Reverse_Image_BW = ~Image_BW_medfilt; figure,imshow(Reverse_Image_BW),title('图象取反'); Image_BW_medfilt2= medfilt2(Reverse_Image_BW,[20 20]); figure,imshow(Image_BW_medfilt2),title('第二次中值滤波的二值化图像'); [Label, Number]=bwlabel(Image_BW_medfilt,8);Number [Label, Number]=bwlabel(Image_BW_medfilt2,8);Number
随着不确定性研究的深入，越来越多的科学家相信，不确定性是这个世界的魅力所在，只有不确定性本身才是确定的。在众多的不确定性中，随机性和模糊性是最基本的。针对概率论和模糊数学在处理不确定性方面的不足，1995年我国工程院院士李德毅教授在概率论和模糊数学的基础上提出了云的概念，并研究了模糊性和随机性及两者之间的关联性。自李德毅院士等人提出云模型至今，云模型已成功的应用到自然语言处理、数据挖掘、 设是一个普通集合。 , 称为论域。关于论域中的模糊集合，是指对于任意元素都存在一个有稳定倾向的随机数，叫做对的隶属度。如果论域中的元素是简单有序的，则可以看作是基础变量，隶属度在上的分布叫做隶属云；如果论域中的元素不是简单有序的，而根据某个法则，可将映射到另一个有序的论域上，中的一个且只有一个和对应，则为基础变量，隶属度在上的分布叫做隶属云[1] 。 数字特征 云模型表示自然语言中的基元——语言值，用云的数字特征
——期望Ex，熵En和超熵He表示语言值的数学性质[3] 。 期望 Ex：云滴在论域空间分布的期望，是最能够代表定性概念的点，是这个概念量化的最典型样本。 熵 En：“熵”这一概念最初是作为描述热力学的一个状态参量，此后又被引入统计物理学、信息论、复杂系统等，用以度量不确定的程度。在云模型中，熵代表定性概念的可度量粒度，熵越大，通常概念越宏观，也是定性概念不确定性的度量，由概念的随机性和模糊性共同决定。一方面, En是定性概念随机性的度量，反映了能够代表这个定性概念的云滴的离散程度；另一方面，又是定性概念亦此亦彼性的度量，反映了在论域空间可被概念接受的云滴的取值范围。用同一个数字特征来反映随机性和模糊性，也必然反映他们之间的关联性。 超熵 He：熵的不确定性度量，即熵的熵，由熵的随机性和模糊性共同决定。反映了每个数值隶属这个语言值程度的凝聚性，即云滴的凝聚程度。超熵越大，云的离散程度越大，隶属度的随机性也随之增大，云的厚度也越大。 1.绘制云图 Ex=18
正弦波的源程序： (一)，用到的函数 1,f2t函数 function x=f2t(X) global dt df t f T N %x=f2t(X) %x为时域的取样值矢量 %X为x的傅氏变换 %X与x长度相同并为2的整幂 %本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)]; x=ifft(X)/dt; end 2,t2f函数。 function X=t2f(x) global dt df N t f T %X=t2f(x) %x为时域的取样值矢量 %X为x的傅氏变换 %X与x长度相同,并为2的整幂。 %本函数需要一个全局变量dt(时域取样间隔) H=fft(x); X=[H(N/2+1:N),H(1:N/2)]*dt; end (二),主程序。 1，%（1）绘出正弦信号波形及频谱 global dt df t f N close all k=input('取样点数＝2^k, k取10左右'); if isempty(k), k=10; end f0=input('f0=取1(kz)左右'); if isempty(f0), f0=1; end N=2^k; dt=0.01; %ms df=1/(N*dt); %KHz T=N*dt; %截短时间
Bs=N*df/2; %系统带宽 f=[-Bs+df/2:df:Bs]; %频域横坐标 t=[-T/2+dt/2:dt:T/2]; %时域横坐标 s=sin(2*pi*f0*t); %输入的正弦信号 S=t2f(s); %S是s的傅氏变换 a=f2t(S); %a是S的傅氏反变换 a=real(a); as=abs(S); subplot(2,1,1) %输出的频谱 plot(f,as,'b'); grid axis([-2*f0,+2*f0,min(as),max(as)]) xlabel('f (KHz)') ylabel('|S(f)| (V/KHz)') %figure(2) subplot(2,1,2) plot(t,a,'black') %输出信号波形画图grid axis([-2/f0,+2/f0,-1.5,1.5]) xlabel('t(ms)') ylabel('a(t)(V)') gtext('频谱图') 最佳基带系统的源程序： (一)，用到的函数 f2t函数和t2f函数。代码>> (二),主程序 globaldt t f df N T close all clear Eb_N0 Pe k=input('取样点数＝2^k, k取13左右'); if isempty(k), k=13; end z=input('每个信号取样点数＝2^z, z
相关主题
文本预览