偏微分方程数值解法的MATLAB源码
[原创]偏微分方程数值解法的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;
matlab源代码实例
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函数源代码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
通信原理基于matlab的计算机仿真_源代码
例错误！文档中没有指定样式的文字。-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
matlab语音识别系统(源代码)最新版
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语者识别的概念 语者识别就是根据说话人的语音信号来判别说话人的身份。语音是人的自然属性之一，由于说话人发音器官的生理差异以及后天形成的行为差异，每个人的语音都带有强烈的个人色彩，这就使得通过分析语音信号来识别说话人成为可能。用语音来鉴别说话人的身份有着许多独特的优点，如语音是人的固有的特征，不会丢失或遗忘；语音信号的采集方便，系统设备成本低；利用电话网络还可实现远程客户服务等。因此，近几年来，说话人识别越来越多的受到人们的重视。与其他生物识别技术如指纹识别、手形识别等相比较，说话人识别不仅使用方便，而且属于非接触性，容易被用户接受，并且在已有的各种生物特征识别技术中，是唯一可以用作远程验证的识别技术。因此，说话人识别的应用前景非常广泛：今天，说话人识别技术已经关系到多学科的研究领域，不同领域中的进步都对说话人识别的发展做出了贡献。说话人识别技术是集声学、语言学、计算机、信息处理和人工智能等诸多领域的一项综合技术，应用需求将十分广阔。在吃力语音信号的时候如何提取信号中关键的成分尤为重要。语音信号的特征参数的好坏直接导致了辨别的准确性。
数字图像处理及matlab实现源代码【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);
偏微分方程数值解法试题与答案
一．填空（1553=?分） 1．若步长趋于零时，差分方程的截断误差0→lm R ，则差分方程的解lm U 趋近于微分方 程的解lm u . 此结论_______（错或对）； 2．一阶Sobolev 空间{} )(,,),()(21 Ω∈''=ΩL f f f y x f H y x 关于内积=1),( g f _____________________是Hilbert 空间； 3．对非线性（变系数）差分格式，常用 _______系数法讨论差分格式的_______稳定性； 4．写出3 x y =在区间]2,1[上的两个一阶广义导数：_________________________________， ________________________________________； 5．隐式差分格式关于初值是无条件稳定的. 此结论_______（错或对）。 二．（13分）设有椭圆型方程边值问题 用1.0=h 作正方形网格剖分 。 （1）用五点菱形差分格式将微分方程在内点离散化； （2）用截断误差为)(2 h O 的差分法将第三边界条件离散化； （3）整理后的差分方程组为 三．（12）给定初值问题 x u t u ??=?? ， ()10,+=x x u 取时间步长1.0=τ，空间步长2.0=h 。试合理选用一阶偏心差分格式（最简显格式）, 并以此格式求出解函数),(t x u 在2.0,2.0=-=t x 处的近似值。 1．所选用的差分格式是： 2．计算所求近似值： 四．（12分）试讨论差分方程 ()h a h a r u u r u u k l k l k l k l ττ + - = -+=++++11,111 1 逼近微分方程 0=??+??x u a t u 的截断误差阶R 。 思路一：将r 带入到原式，展开后可得格式是在点（l+1/2,k+1/2）展开的。 思路二：差分格式的用到的四个点刚好是矩形区域的四个顶点，可由此构造中心点的差分格 式。
基于MATLAB的潮流计算源程序代码(优.选)
%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算********** 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;
数字图像处理_旋转与幅度谱(含MATLAB代码)
数字图像处理实验一 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)
偏微分方程数值解期末试题及答案(内容参考)
偏微分方程数值解试题(06B) 参考答案与评分标准 信息与计算科学专业 一（10分）、设矩阵A 对称，定义)(),(),(2 1 )(n R x x b x Ax x J ∈-= ，)()(0x x J λλ?+=.若0)0('=?,则称称0x 是)(x J 的驻点（或稳定点）.矩阵A 对称（不必正定），求证0x 是)(x J 的驻点的充要条件是：0x 是方程组 b Ax =的解 解: 设n R x ∈0是)(x J 的驻点,对于任意的n R x ∈,令 ),(2 ),()()()(2 000x Ax x b Ax x J x x J λλλλ?+ -+=+=, (3分) 0)0('=?,即对于任意的n R x ∈,0),(0=-x b Ax ,特别取b Ax x -=0,则有 0||||),(2000=-=--b Ax b Ax b Ax ,得到b Ax =0. (3分) 反之,若 n R x ∈0满足 b Ax =0,则对于任意的 x ,)(),(2 1 )0()1()(00x J x Ax x x J >+ ==+??,因此0x 是)(x J 的最小值点. (4分) 评分标准:)(λ?的展开式3分, 每问3分,推理逻辑性1分 二（10分）、 对于两点边值问题：????? ==∈=+-=0 )(,0)() ,()(' b u a u b a x f qu dx du p dx d Lu 其中]),([,0]),,([,0)(min )(]),,([0min ] ,[1b a H f q b a C q p x p x p b a C p b a x ∈≥∈>=≥∈∈ 建立与上述两点边值问题等价的变分问题的两种形式：求泛函极小的Ritz 形式和 Galerkin 形式的变分方程。 解: 设}0)(),,(|{11 =∈=a u b a H u u H E 为求解函数空间,检验函数空间.取),(1 b a H v E ∈,乘方程两端,积分应用分部积分得到 (3分) )().(),(v f fvdx dx quv dx dv dx du p v u a b a b a ==+=??,),(1 b a H v E ∈? 即变分问题的Galerkin 形式. (3分)
matlab图像处理代码
附录 ＭＡＴＬＡＢ图像处理命令 １．ａｐｐｌｙｌｕｔ 功能： 在二进制图像中利用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);
偏微分方程数值解法
一、 问题 用有限元方法求下面方程的数值解 2 u u u f t ?-?+=? in (]0,T Ω? 0u = on []0,T ?Ω? ()00,u x u = in Ω 二、 问题分析 第一步 利用Green 公式，求出方程的变分形式 变分形式为：求()()21 00,;u L T H ∈Ω，使得 ()())(2 ,,,,u v u v u v f v t ???+??+= ???? ()10v H ?∈Ω （*） 以及 ()00,u x u =. 第二步 对空间进行离散，得出半离散格式 对区域Ω进行剖分，构造节点基函数，得出有限元子空间：()12,,,h NG V span ???=???，则（*）的Galerkin 逼近为： []0,t T ?∈，求()()1 0,h h u t x V H ∈?Ω，使得 ()()()()() () )(2 ,,,,h h h h h h h d u t v u t v u t v f v dt +??+= h h v V ?∈ （**） 以及()0,0h h u u =，0,h u 为初始条件0u 在h V 中的逼近，设0,h u 为0u 在h V 中的插值. 则0t ?≥，有()()1 N G h i i i u t t ξ? == ∑，0,h u =01 N G i i i ξ?=∑，代人（**）即可得到一常微分方程组. 第三步 进一步对时间进行离散，得到全离散的逼近格式 对 du dt 用差分格式.为此把[]0,T 等分为n 个小区间[]1,i i t t -，其长度1i i T t t t n -?=-= ,n t T =. 这样把求i t 时刻的近似记为i h u ，0 h u 是0u 的近似.这里对（**）采用向后的欧拉格式，即 ()()() () )(2 11 11 1 ,,,,i i i i h h h h h h h i h u u v u v u v f v t ++++-+??+ = ? h h v V ?∈ (***) i=0,1,2…,n-1. 0 h u =0,h u 由于向后欧拉格式为隐式格式且含有非线性项，故相邻两时间步之间采用牛顿迭代，即：
最常用的matlab图像处理的源代码
最常用的一些图像处理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);
数字图像处理matlab代码
一、编写程序完成不同滤波器的图像频域降噪和边缘增强的算法并进行比较，得出结论。 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赋初值
基本粒子群算法的matlab源程序
主函数源程序（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)
BP神经网络matlab源程序代码
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');
图像处理实例(含Matlab代码)
信号与系统实验报告——图像处理 学院：信息科学与工程学院 专业：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
Matlab源程序代码
正弦波的源程序： (一)，用到的函数 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 图像处理matlab源程序 1.图像反转 MATLAB程序实现如下： I=imread('xian.bmp'); J=double(I); J=-J+(256-1); %图像反转线性变换 H=uint8(J); subplot(1,2,1),imshow(I); subplot(1,2,2),imshow(H); 2.灰度线性变换 MATLAB程序实现如下： I=imread('xian.bmp'); subplot(2,2,1),imshow(I); title('原始图像'); axis([50,250,50,200]); axis on; %显示坐标系 I1=rgb2gray(I); subplot(2,2,2),imshow(I1); title('灰度图像'); axis([50,250,50,200]); axis on; %显示坐标系 J=imadjust(I1,[0.1 0.5],[]); %局部拉伸，把[0.1 0.5]内的灰度拉伸为[0 1] subplot(2,2,3),imshow(J); title('线性变换图像[0.1 0.5]'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 K=imadjust(I1,[0.3 0.7],[]); %局部拉伸，把[0.3 0.7]内的灰度拉伸为[0 1] subplot(2,2,4),imshow(K); title('线性变换图像[0.3 0.7]'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 3.非线性变换 MATLAB程序实现如下： I=imread('xian.bmp'); I1=rgb2gray(I); subplot(1,2,1),imshow(I1); title('灰度图像'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 J=double(I1);图像处理matlab源程序