灰色预测matlab程序

2014.4.11~4.20 预测价格为 5.9535

>>%灰色预测模型键搏程序

clear

syms a b

c=[a b]'

A=[5.76 5.92 5.91 5.92 5.93 5.93 5.93 5.94 5.94 5.95 5.95 5.94 ]% 原始序列

B=cumsum(A)%累加

n=length(A)

for i=1:(n-1)

C(i)=(B(i)+B(i+1))/2

end

%计算待定参数

D=A

D(1)=[]

D=D'

E=[-Cones(1,n-1)]

c=inv(E*E')*E*D

c=c'

a=c(1)

b=c(2)

%预测 往后预测1个数据（搏亮胡即4.11-4.20价格）

F=[]F(1)=A(1)

for i=2:(n+1)

F(i)=(A(1)-b/a)/exp(a*(i-1))+b/a

end

G=[]G(1)=A(1)

for i=2:(n+1)

G(i)=F(i)-F(i-1)

end

t1=2002:2013

t2=2002:2014

G

plot(t1,A,'基拦o',t2,G)

G =

5.76005.91555.91895.92245.92585.92935.93275.93625.93965.94315.94655.95005.9535

Matlab上有CANNY算子的库函数啊，直接调用就行了。

我这有VC++的边缘检测算法，很长的。稍微改一下就可以缺宴森用在Matlab上。

/ 一维高斯祥中分布函数，用于平滑函数中生成的高斯滤波系数

void CFunction::CreatGauss(double sigma, double **pdKernel, int *pnWidowSize)

{

LONG i

//数组中心点

int nCenter

//数组中一点到中心点距离

double dDis

//中间变量

double dValue

double dSum

dSum = 0

// [-3*sigma,3*sigma] 以内数据，会覆盖绝大部分滤波系数

*pnWidowSize = 1+ 2*ceil(3*sigma)

nCenter = (*pnWidowSize)/2

*pdKernel = new double[*pnWidowSize]

//生伏亩成高斯数据

for(i=0i<(*pnWidowSize)i++)

{

dDis = (double)(i - nCenter)

dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma))/(sqrt(2*3.1415926)*sigma)

(*pdKernel)[i] = dValue

dSum+=dValue

}

//归一化

for(i=0i<(*pnWidowSize)i++)

{

(*pdKernel)[i]/=dSum

}

}

//用高斯滤波器平滑原图像

void CFunction::GaussianSmooth(SIZE sz, LPBYTE pGray, LPBYTE pResult, double sigma)

{

LONG x, y

LONG i

//高斯滤波器长度

int nWindowSize

//窗口长度

int nLen

//一维高斯滤波器

double *pdKernel

//高斯系数与图像数据的点乘

double dDotMul

//滤波系数总和

double dWeightSum

double *pdTemp

pdTemp = new double[sz.cx*sz.cy]

//产生一维高斯数据

CreatGauss(sigma, &pdKernel, &nWindowSize)

nLen = nWindowSize/2

//x方向滤波

for(y=0y<sz.cyy++)

{

for(x=0x<sz.cxx++)

{

dDotMul = 0

dWeightSum = 0

for(i=(-nLen)i<=nLeni++)

{

//判断是否在图像内部

if((i+x)>=0 &&(i+x)<sz.cx)

{

dDotMul+=(double)pGray[y*sz.cx+(i+x)] * pdKernel[nLen+i]

dWeightSum += pdKernel[nLen+i]

}

}

pdTemp[y*sz.cx+x] = dDotMul/dWeightSum

}

}

//y方向滤波

for(x=0x<sz.cxx++)

{

for(y=0y<sz.cyy++)

{

dDotMul = 0

dWeightSum = 0

for(i=(-nLen)i<=nLeni++)

{

if((i+y)>=0 &&(i+y)<sz.cy)

{

dDotMul += (double)pdTemp[(y+i)*sz.cx+x]*pdKernel[nLen+i]

dWeightSum += pdKernel[nLen+i]

}

}

pResult[y*sz.cx+x] = (unsigned char)(int)dDotMul/dWeightSum

}

}

delete []pdKernel

pdKernel = NULL

delete []pdTemp

pdTemp = NULL

}

// 方向导数,求梯度

void CFunction::Grad(SIZE sz, LPBYTE pGray,int *pGradX, int *pGradY, int *pMag)

{

LONG y,x

//x方向的方向导数

for(y=1y<sz.cy-1y++)

{

for(x=1x<sz.cx-1x++)

{

pGradX[y*sz.cx +x] = (int)( pGray[y*sz.cx+x+1]-pGray[y*sz.cx+ x-1] )

}

}

//y方向方向导数

for(x=1x<sz.cx-1x++)

{

for(y=1y<sz.cy-1y++)

{

pGradY[y*sz.cx +x] = (int)(pGray[(y+1)*sz.cx +x] - pGray[(y-1)*sz.cx +x])

}

}

//求梯度

//中间变量

double dSqt1

double dSqt2

for(y=0y<sz.cyy++)

{

for(x=0x<sz.cxx++)

{ //二阶范数求梯度

dSqt1 = pGradX[y*sz.cx + x]*pGradX[y*sz.cx + x]

dSqt2 = pGradY[y*sz.cx + x]*pGradY[y*sz.cx + x]

pMag[y*sz.cx+x] = (int)(sqrt(dSqt1+dSqt2)+0.5)

}

}

}

//非最大抑制

void CFunction::NonmaxSuppress(int *pMag, int *pGradX, int *pGradY, SIZE sz, LPBYTE pNSRst)

{

LONG y,x

int nPos

//梯度分量

int gx

int gy

//中间变量

int g1,g2,g3,g4

double weight

double dTmp,dTmp1,dTmp2

//设置图像边缘为不可能的分界点

for(x=0x<sz.cxx++)

{

pNSRst[x] = 0

//pNSRst[(sz.cy-1)*sz.cx+x] = 0

pNSRst[sz.cy-1+x] = 0

}

for(y=0y<sz.cyy++)

{

pNSRst[y*sz.cx] = 0

pNSRst[y*sz.cx + sz.cx-1] = 0

}

for(y=1y<sz.cy-1y++)

{

for(x=1x<sz.cx-1x++)

{ //当前点

nPos = y*sz.cx + x

//如果当前像素梯度幅度为0，则不是边界点

if(pMag[nPos] == 0)

{

pNSRst[nPos] = 0

}

else

{ //当前点的梯度幅度

dTmp = pMag[nPos]

//x,y方向导数

gx = pGradX[nPos]

gy = pGradY[nPos]

//如果方向导数y分量比x分量大，说明导数方向趋向于y分量

if(abs(gy) >abs(gx))

{

//计算插值比例

weight = fabs(gx)/fabs(gy)

g2 = pMag[nPos-sz.cx]

g4 = pMag[nPos+sz.cx]

//如果x,y两个方向导数的符号相同

//C 为当前像素，与g1-g4 的位置关系为：

//g1 g2

// C

// g4 g3

if(gx*gy>0)

{

g1 = pMag[nPos-sz.cx-1]

g3 = pMag[nPos+sz.cx+1]

}

//如果x,y两个方向的方向导数方向相反

//C是当前像素，与g1-g4的关系为：

// g2 g1

// C

// g3 g4

else

{

g1 = pMag[nPos-sz.cx+1]

g3 = pMag[nPos+sz.cx-1]

}

}

//如果方向导数x分量比y分量大，说明导数的方向趋向于x分量

else

{

//插值比例

weight = fabs(gy)/fabs(gx)

g2 = pMag[nPos+1]

g4 = pMag[nPos-1]

//如果x,y两个方向的方向导数符号相同

//当前像素C与 g1-g4的关系为

// g3

// g4 C g2

// g1

if(gx * gy >0)

{

g1 = pMag[nPos+sz.cx+1]

g3 = pMag[nPos-sz.cx-1]

}

//如果x,y两个方向导数的方向相反

// C与g1-g4的关系为

// g1

// g4 C g2

// g3

else

{

g1 = pMag[nPos-sz.cx+1]

g3 = pMag[nPos+sz.cx-1]

}

}

//利用 g1-g4 对梯度进行插值

{

dTmp1 = weight*g1 + (1-weight)*g2

dTmp2 = weight*g3 + (1-weight)*g4

//当前像素的梯度是局部的最大值

//该点可能是边界点

if(dTmp>=dTmp1 &&dTmp>=dTmp2)

{

pNSRst[nPos] = 128

}

else

{

//不可能是边界点

pNSRst[nPos] = 0

}

}

}

}

}

}

// 统计pMag的直方图，判定阈值

void CFunction::EstimateThreshold(int *pMag, SIZE sz, int *pThrHigh, int *pThrLow, LPBYTE pGray,

double dRatHigh, double dRatLow)

{

LONG y,x,k

//该数组的大小和梯度值的范围有关，如果采用本程序的算法

//那么梯度的范围不会超过pow(2,10)

int nHist[1024]

//可能边界数

int nEdgeNum

//最大梯度数

int nMaxMag

int nHighCount

nMaxMag = 0

//初始化

for(k=0k<1024k++)

{

nHist[k] = 0

}

//统计直方图,利用直方图计算阈值

for(y=0y<sz.cyy++)

{

for(x=0x<sz.cxx++)

{

if(pGray[y*sz.cx+x]==128)

{

nHist[pMag[y*sz.cx+x]]++

}

}

}

nEdgeNum = nHist[0]

nMaxMag = 0

//统计经过“非最大值抑制”后有多少像素

for(k=1k<1024k++)

{

if(nHist[k] != 0)

{

nMaxMag = k

}

//梯度为0的点是不可能为边界点的

//经过non-maximum suppression后有多少像素

nEdgeNum += nHist[k]

}

//梯度比高阈值*pThrHigh 小的像素点总书目

nHighCount = (int)(dRatHigh * nEdgeNum + 0.5)

k=1

nEdgeNum = nHist[1]

//计算高阈值

while((k<(nMaxMag-1)) &&(nEdgeNum <nHighCount))

{

k++

nEdgeNum += nHist[k]

}

*pThrHigh = k

//低阈值

*pThrLow = (int)((*pThrHigh) * dRatLow + 0.5)

}

//利用函数寻找边界起点

void CFunction::Hysteresis(int *pMag, SIZE sz, double dRatLow, double dRatHigh, LPBYTE pResult)

{

LONG y,x

int nThrHigh,nThrLow

int nPos

//估计TraceEdge 函数需要的低阈值，以及Hysteresis函数使用的高阈值

EstimateThreshold(pMag, sz,&nThrHigh,&nThrLow,pResult,dRatHigh,dRatLow)

//寻找大于dThrHigh的点，这些点用来当作边界点，

//然后用TraceEdge函数跟踪该点对应的边界

for(y=0y<sz.cyy++)

{

for(x=0x<sz.cxx++)

{

nPos = y*sz.cx + x

//如果该像素是可能的边界点，并且梯度大于高阈值，

//该像素作为一个边界的起点

if((pResult[nPos]==128) &&(pMag[nPos] >= nThrHigh))

{

//设置该点为边界点

pResult[nPos] = 255

TraceEdge(y,x,nThrLow,pResult,pMag,sz)

}

}

}

//其他点已经不可能为边界点

for(y=0y<sz.cyy++)

{

for(x=0x<sz.cxx++)

{

nPos = y*sz.cx + x

if(pResult[nPos] != 255)

{

pResult[nPos] = 0

}

}

}

}

//根据Hysteresis 执行的结果，从一个像素点开始搜索，搜索以该像素点为边界起点的一条边界的

//一条边界的所有边界点，函数采用了递归算法

// 从（x,y)坐标出发，进行边界点的跟踪，跟踪只考虑pResult中没有处理并且可能是边界

// 点的像素（=128），像素值为0表明该点不可能是边界点，像素值为255表明该点已经是边界点

void CFunction::TraceEdge(int y, int x, int nThrLow, LPBYTE pResult, int *pMag, SIZE sz)

{

//对8邻域像素进行查询

int xNum[8] = {1,1,0,-1,-1,-1,0,1}

int yNum[8] = {0,1,1,1,0,-1,-1,-1}

LONG yy,xx,k//循环变量

for(k=0k<8k++)

{

yy = y+yNum[k]

xx = x+xNum[k]

if(pResult[640 * (479 - yy)+xx]==128 &&pMag[640 * (479 - yy)+xx]>=nThrLow )

{

//该点设为边界点

pResult[640 * (479 - yy)+xx] = 255

//以该点为中心再进行跟踪

TraceEdge(yy,xx,nThrLow,pResult,pMag,sz)

}

}

}

// Canny算子

BOOL CFunction::Canny(LPBYTE m_pDibData,CPoint ptLeft, CPoint ptRight , double sigma, double dRatLow, double dRatHigh)

{

BYTE* m_Newdata//每一步处理后的图像数据

m_Newdata = (BYTE*)malloc(maxImage)

memcpy(m_Newdata,(BYTE *)m_pDibData,maxImage)

//经过抑制局部像素非最大值的处理后的数据

BYTE* pResult//每一步处理后的图像数据

pResult = (BYTE*)malloc(maxImage)

memcpy(pResult,(BYTE *)m_pDibData,maxImage)

int pointy,pointx,m,n,i=0

long Position

int GradHori

int GradVert

//存储结构元素的数组

BYTE array[9]={0}

//设定两个阈值

int nThrHigh,nThrLow

//梯度分量

int gx

int gy

//中间变量

int g1,g2,g3,g4

double weight

double dTmp,dTmp1,dTmp2

int Width,Higth

Width=ptRight.x-ptLeft.x+1

Higth=ptRight.y-ptLeft.y+1

CSize sz=CSize(Width,Higth)

//x方向导数的指针

int *pGradX= new int[maxImage]

memset(pGradX,0,maxImage)

//y方向

int *pGradY

pGradY = new int [maxImage]

memset(pGradY,0,maxImage)

//梯度的幅度

int *pGradMag

pGradMag = new int [maxImage]

//对pGradMag进行初始化

for (pointy = 0pointy <480pointy++)

{

for (pointx = 0pointx <640 pointx++)

{

Position=640 * (479 - pointy)+pointx

pGradMag[Position]=m_pDibData[Position]

}

}

//第一步进行高斯平滑器滤波

//进入循环,使用3*3的结构元素，处理除去第一行和最后一行以及第一列和最后一列。

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{

Position=640 * (479 - pointy)+pointx

for (m = 0m <3m++)

{

for (n = 0n <3n++)

{

array[m*3+n]=m_pDibData[Position+640*(1-m)+n-1]

}

}

GradHori=abs(array[0]+2*array[1]+array[2]+2*array[3]+4*array[4]+2*array[5]+array[6]+2*array[7]+array[8])

GradHori=(int)(0.0625*GradHori+0.5)

if (GradHori>255)

{

m_Newdata[Position]=255

}

else

m_Newdata[Position]=GradHori

}

}

//第二步用一阶偏导的有限差分来计算梯度的幅值和方向

//x方向的方向导数

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{

pGradX[pointy*Width +pointx]=(int)(m_Newdata[pointy*Width +pointx+1]- m_Newdata[pointy*Width +pointx-1] )

}

}

//y方向方向导数

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

pGradY[pointy*Width +pointx] = (int)(m_Newdata[(pointy+1)*Width +pointx] - m_Newdata[(pointy-1)*Width +pointx])

}

}

//求梯度

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{

Position=640 * (479 - pointy)+pointx

for (m = 0m <3m++)

{

for (n = 0n <3n++)

{

array[m*3+n]=m_Newdata[Position+640*(1-m)+n-1]

}

}

GradHori=abs((-1)*array[0]+(-2)*array[3]+2*array[7]+array[8])

GradVert=abs((-1)*array[0]-2*array[1]+2*array[5]+array[8])

GradHori =(int)((float)sqrt(pow(GradHori,2)+pow(GradVert,2))+0.5)

pGradMag[Position]=GradHori

}

}

//针对第一行的像素点及最后一行的像素点

for (pointx = ptLeft.xpointx <= ptRight.xpointx++)

{

Position=640 * (479 - ptLeft.y)+pointx

pGradMag[Position]=0

Position=640 * (479 - ptRight.y)+pointx

pGradMag[Position]=0

}

//针对第一列以及最后一列的像素点

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

Position=640 * (479 - pointy)+ptLeft.x

pGradMag[Position]=0

Position=640 * (479 - pointy)+ptRight.x

pGradMag[Position]=0

}

//第三步进行抑制梯度图中的非局部极值点的像素

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{ //当前点

Position=640 * (479 - pointy)+pointx

//如果当前像素梯度幅度为0，则不是边界点

if(pGradMag[Position] == 0)

{

pGradMag[Position] = 0

}

else

{ //当前点的梯度幅度

dTmp = pGradMag[Position]

//x,y方向导数

gx = pGradX[Position]

gy = pGradY[Position]

//如果方向导数y分量比x分量大，说明导数方向趋向于y分量

if(abs(gy) >abs(gx))

{

//计算插值比例

weight = fabs(gx)/fabs(gy)

g2 = pGradMag[Position-640]

g4 = pGradMag[Position+640]

//如果x,y两个方向导数的符号相同

//C 为当前像素，与g1-g4 的位置关系为：

//g1 g2

//C

//g4 g3

if(gx*gy>0)

{

g1 = pGradMag[Position-640-1]

g3 = pGradMag[Position+640+1]

}

//如果x,y两个方向的方向导数方向相反

//C是当前像素，与g1-g4的关系为：

// g2 g1

// C

// g3 g4

else

{

g1 = pGradMag[Position-640+1]

g3 = pGradMag[Position+640-1]

}

}

//如果方向导数x分量比y分量大，说明导数的方向趋向于x分量

else

{

//插值比例

weight = fabs(gy)/fabs(gx)

g2 = pGradMag[Position+1]

g4 = pGradMag[Position-1]

//如果x,y两个方向的方向导数符号相同

//当前像素C与 g1-g4的关系为

// g3

// g4 C g2

// g1

if(gx * gy >0)

{

g1 = pGradMag[Position+640+1]

g3 = pGradMag[Position-640-1]

}

//如果x,y两个方向导数的方向相反

// C与g1-g4的关系为

// g1

// g4 C g2

// g3

else

{

g1 =pGradMag[Position-640+1]

g3 =pGradMag[Position+640-1]

}

}

//利用 g1-g4 对梯度进行插值

{

dTmp1 = weight*g1 + (1-weight)*g2

dTmp2 = weight*g3 + (1-weight)*g4

//当前像素的梯度是局部的最大值

//该点可能是边界点

if(dTmp>=dTmp1 &&dTmp>=dTmp2)

{

pResult[Position] = 128

}

else

{

//不可能是边界点

pResult[Position] = 0

}

}

}

}

}

//第四步根据梯度计算及经过非最大值得印制后的结果设定阈值

//估计TraceEdge 函数需要的低阈值，函数使用的高阈值

EstimateThreshold(pGradMag, sz,&nThrHigh,&nThrLow,pResult,dRatHigh,dRatLow)

//寻找大于dThrHigh的点，这些点用来当作边界点，

//然后用TraceEdge函数跟踪该点对应的边界

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{

Position=640 * (479 - pointy)+pointx

//如果该像素是可能的边界点，并且梯度大于高阈值，

//该像素作为一个边界的起点

if((pResult[Position]==128) &&(pGradMag[Position] >= nThrHigh))

{

//设置该点为边界点

pResult[Position] = 255

TraceEdge(pointy,pointx,nThrLow,pResult,pGradMag,sz)

}

}

}

//其他点已经不可能为边界点

for (pointy = ptLeft.y+1pointy <= ptRight.y-1pointy++)

{

for (pointx = ptLeft.x+1pointx <= ptRight.x-1pointx++)

{

Position=640 * (479 - pointy)+pointx

if(pResult[Position] != 255)

{

pResult[Position] = 0

}

}

}

//计算方向导数和梯度的幅度

// Grad(sz,pGaussSmooth,pGradX,pGradY,pGradMag)

//应用非最大抑制

// NonmaxSuppress(pGradMag,pGradX,pGradY,sz,pResult)

//应用Hysteresis，找到所有边界

// Hysteresis(pGradMag,sz,dRatLow,dRatHigh,pResult)

memcpy(m_pDibData,(BYTE *)pResult,maxImage)

delete[] pResult

pResult = NULL

delete[] pGradX

pGradX = NULL

delete[] pGradY

pGradY = NULL

delete[] pGradMag

pGradMag = NULL

delete[] m_Newdata

m_Newdata = NULL

return true

}

1、打开MATLAB软件，准备开始作图。

方法一：

1、在命令窗口物蚂，编写离散点向量，例如离散点(X，Y)，X= [1 3 5 7 2 9 3 6 2 8]，Y= [2 1 4 5 6 15 2 5 7 10]，如果我们直接采用Plot作图，画出来的是折线。

2、对Plot显示样式进行设置，如：plot(X,Y,'k*')，'k'表示显示黑色颜色，'*'表示点为星号显示。

方法二：

1、也比较方便，清除之前在命令行的代码，直接在命令行输入：clc，然后按回车键。

2、重新在命令行输入离散点，采用函数scatter进行画离散点，如：scatter(X,Y,'r')，'r'表示离散点显示为红色。

3、比较一下两种方法是否一样，用hold on ,使两个图在一副图上面显示出来。发现标记的离散点是在同一个位置的！

在使锋伍用Python绘制图表前，我们需要先安装两个库文件numpy和matplotlib。

Numpy是Python开源的数值计算扩展，可用来存储和罩基埋处理大型矩阵，比Python自身数据结构要高效；matplotlib是一个Python的图像框架，使用其绘制出来的图形效果和MATLAB下绘制的图形类似。

《精通matlab数字图像处理与识别》将理论知识、科学研究和工程实践有机结合起来，内容涉及数字图像处理和识别技术的方方面面，包括图像的点运算、几何变换、空域和频域滤波、小波变换、图像复原、形态学处理、图像分割以及图像特征提取的相关内容；同时对于机器视觉进行了前导性的探究，重点介绍了两种目前在工程技术领域非常流行的分类技术——人工神经网络(ann)和支持向量机(svm)，并在人脸识别这样的热点问题中结束本书。

《精通matlab数字图像处理与识别》结构紧凑，内容深入浅出，讲解图文并茂，适合计算机、通信和自动化等相关专业的本科生、研究生，以及工作在图像处理和识别领域一线的广大工程技术人员参考使用。

---