问题描述：
及最后预测过程
线性预测出的波形幅度是不是会衰减
问题解答：
Syntax
r = rlevinson(a,efinal)
[r,u] = rlevinson(a,efinal)
[r,u,k] = rlevinson(a,efinal)
[r,u,k,e] = rlevinson(a,efinal)
Description
The reverse Levinson-Durbin recursion implements the step-down algorithm for solving the following symmetric Toeplitz system of linear equations for r, where r = [r(1) L r(p+1)] and r(i)* denotes the complex conjugate of r(i).
r = rlevinson(a,efinal) solves the above system of equations for r given vector a, where a = [1 a(2) L a(p+1)]. In linear prediction applications, r represents the autocorrelation sequence of the input to the prediction error filter, where r(1) is the zero-lag element. The figure below shows the typical filter of this type, where H(z) is the optimal linear predictor, x(n) is the input signal, is the predicted signal, and e(n) is the prediction error.
Input vector a represents the polynomial coefficients of this prediction error filter in descending powers of z.
The filter must be minimum phase to generate a valid autocorrelation sequence. efinal is the scalar prediction error power, which is equal to the variance of the prediction error signal, σ2(e).
[r,u] = rlevinson(a,efinal) returns upper triangular matrix U from the UDU* decomposition
where
and E is a diagonal matrix with elements returned in output e (see below). This decomposition permits the efficient evaluation of the inverse of the autocorrelation matrix, R-1.
Output matrix u contains the prediction filter polynomial, a, from each iteration of the reverse Levinson-Durbin recursion
where ai(j) is the jth coefficient of the ith order prediction filter polynomial (i.e., step i in the recursion). For example, the 5th order prediction filter polynomial is
a5 = u(5:-1:1,5)'
Note that u(p+1:-1:1,p+1)' is the input polynomial coefficient vector a.
[r,u,k] = rlevinson(a,efinal) returns a vector k of length (p+1) containing the reflection coefficients. The reflection coefficients are the conjugates of the values in the first row of u.
k = conj(u(1,2:end))
[r,u,k,e] = rlevinson(a,efinal) returns a vector of length p+1 containing the prediction errors from each iteration of the reverse Levinson-Durbin recursion: e(1) is the prediction error from the first-order model, e(2) is the prediction error from the second-order model, and so on.
These prediction error values form the diagonal of the matrix E in the UDU* decomposition of R-1.
#include <iostream.h>
#include <fstream.h>
#include <math.h>
const int order=8;
const int N=256;
const int RequestN=1024;
const float PI=3.14159;
float signal[N];//输入已测信号
float D[order+1],V[order+1][order+1],A[order+1],Y[order+1];
float cov(int ,int);
void getYi();
void getAi();
float getDi(int );
float getVij(int ,int );
void initData()
{
for (int i=0;i<=order;i++)
{
D[i]=0;
A[i]=0;
Y[i]=0;
for (int j=0;j<=order;j++)
V[i][j]=0;
}
}
float cov(int i,int j)
{
int n=i>j?i:j;
float co=0;
while (n<N)
{
co+=signal[n-i-1]*signal[n-j-1];
n++;
}
return co;
}
float getVij(int i,int j)
{
float sum=0;
for (int k=1;k<j;k++)sum+=V[i][k]*D[k]*V[j][k];
V[i][j]=(cov(i,j)-sum)/D[j];
}
float getDi(int i)
{
float sum=0;
for (int k=1;k<i;k++)
sum+=V[i][k]*V[i][k]*D[k];
return cov(i,i)-sum;
}
void getAi()
{
A[order]=Y[order]/D[order];
int i,j;
for (i=order-1;i>=1;i--)
{
float sum=0;
for (j=i+1;j<=order;j++)
sum+=V[j][i]*A[j];
A[i]=Y[i]/D[i]-sum;
}
}
void getYi()
{
Y[1]=cov(1,0);
for (int i=2;i<=order;i++)
{
float sum=0;
for (int j=1;j<i-1;j++)
sum+=V[i][j]*Y[j];
Y[i]=cov(i,0)-sum;
}
}
int main(void)
{
initData();
for (int i=0;i<N;i++)
signal[i]=10*sin(PI*4/N*i);
D[1]=cov(1,1);
for (int i=1;i<=order;i++)
{
for (int j=i+1;j<=order;j++)
V[j][i]=getVij(i,j);
if (i+1<=order)
D[i+1]=getDi(i+1);
}
getYi();
getAi();
for (int jj=1;jj<order+1;jj++)
cout<<Y[jj]<<' ';
return 0;
}
协方差的计算机公式为
Sxy=(∑(xi-x^)(yi-y^))/(n-1)，i=0,1,2...,n-1，这里^表示平均值
下面就是协方差法的C实现，并且给了一个main函数进行验证：
#include<conio.h>
#include<stdio.h>
void cox(int x[],int y[],int n,double *p)
{
int i;
double averx,avery,sumx=0,sumy=0;
double sum=0;
for(i=0;i<n;i++)
{
sumx+=x[i];
sumy+=y[i];
}
averx=sumx/n;
avery=sumy/n;
for(i=0;i<n;i++)
sum+=(x[i]-averx)*(y[i]-avery);
*p=sum/(n-1);
}
void output(int array[],int n)
{
int i;
for(i=0;i<n;i++)
printf("%4d",array[i]);
}
int main()
{
int x[5]={1,2,3,4,5};
int y[5]={5,4,3,2,1};
double co;
cox(x,y,5,&co);
printf("\nx:\n");
output(x,5);
printf("\ny:\n");
output(y,5);
printf("\nThe cox of x and y is:%lf\n",co);
getch();
return 1;
}
上海旅游网