[DSP编程] 我的5509A之路

clc; clear all;

%% lms

N = 2048;                               % 信号的取样点数

M = 10;                                  % 抽头的个数

iter = N;                               % 迭代次数

t = 1 : N;

s = 5 * sin(2 * pi * 0.01 * t);         % 源信号

x = s + 0.5 * randn(1, N);              % 输入信号

x = x';

d = s';                                 % 期望信号

e = zeros(iter, 1);            

w = zeros(M, iter);   



rho_max = max(eig(x * x'));   % 输入信号相关矩阵的最大特征值

mu = rand(1) * (2 / rho_max)   ;   % 收敛因子 0 < mu < 2/rho



for n = M + 1 : iter -1                

    w( :, n + 1) = w( :, n) + mu * x(n - 1 : -1 : n - M) * conj(e(n));

    y(n + 1) = w( :, n + 1)' * x(n : -1 : n - M + 1);

    e(n + 1) = d(n + 1, : ) - y(n + 1);

end



%% 画输入输出信号以及输出信号和预期输出信号的误差

figure(1);

subplot(2, 1, 1);

plot(t, x);

grid;

title('输入信号');

xlabel('样点数');

ylabel('幅值');



subplot(2,1,2);

plot(t, y);

grid;

title('输出信号');

xlabel('样点数');

ylabel('幅值');



% 绘制自适应滤波器输出信号,预期输出信号和两者的误差

figure(2) 

subplot(3, 1, 1);

plot(t, d);

grid;

title('预期输出');

xlabel('样点数');

ylabel('幅值');



subplot(3, 1, 2);

plot(t, y);

grid;

title('自适应滤波器输出');

xlabel('样点数');

ylabel('幅值');



subplot(3, 1, 3);

plot(t, e);

grid;

title('误差');

xlabel('样点数');

ylabel('幅值');

clc; clear all;

%% 产生512点输入序列

M = 2;                                               % AR 阶数

N = 512;                                             % 样本数

trials = 100;                                        % 试验次数

a1 = -0.975;                                         % AR parameter

a2 = 0.95;

var = 0.0731;                                        % 噪声方差

v = 0 + sqrt(var).* randn(N, 1);                     % 产生v(n)

u0 = [0, 0];                                         % u(0)初始化

a = [1 a1 a2];

b = 1;

zi = filtic(b, a, u0);             % 初始化滤波器 直接II型

[u, zf] = filter(b, a, v, zi);     % 产生观测信号u(n)



%% 用LMS算法计算w1和w2



mu1 = 0.05;

mu2 = 0.005;

w1 = zeros(M, N);                                   % 初始化权向量

w2 = zeros(M, N); 

w1_trials = zeros(M, N, trials);

w2_trials = zeros(M, N, trials);

                       

e1 = zeros(N, 1);                                   % 初始化误差

e2 = zeros(N, 1);

e1_trials = zeros(N, 1, trials);  

e2_trials = zeros(N, 1, trials); 



d1 = zeros(N, 1);                                   % 初始化期望响应

d2 = zeros(N, 1);



wopt = zeros(2, trials);

Jmin = zeros(1, trials);   

sum_eig = zeros(1, trials);



for i = 1 : trials

     v = 0 + sqrt(var).* randn(N, 1);

     zi = filtic(b, a, u0);              

     [u, zf] = filter(b, a, v, zi);      

     

     uu = zeros(M, N - M + 1);                       

     for k = 1 : (N - M + 1)

         uu( :, k) = u(k + M - 1 : -1 : k).';          

     end

     R = uu * uu' / (N - M + 1);                     % u(n)自相关矩阵

     p = sum(uu(:, 1 : N - M) * diag(u(3 : N)), 2) / (N - M);

     [v, d] = eig(R);

     wopt( :, i) = inv(R) * p;

     Jmin(i)   = R(1, 1) - p' * wopt( :, i);         % Jmin

     sum_eig(i) = d(1, 1) + d(2, 2);                 % 特征值之和

   

     for n = 3 : N -1

        w1( :, n + 1) = w1( :, n) + mu1 * u(n - 1 : -1 : n - 2, :) * conj(e1(n));

        d1(n + 1) = w1( :, n + 1)' * u(n : -1 : n - 1, : );

        e1(n + 1) = u(n + 1, : ) - d1(n + 1);

        

        w2( :, n + 1) = w2( :, n) + mu2 * u(n - 1 : -1 : n - 2, :) * conj(e2(n));

        d2(n + 1) = w2( :, n + 1)' * u(n : -1 : n - 1, : );

        e2(n + 1) = u(n + 1, : ) - d2(n + 1);

     end

     w1_trials ( :, :, i) = w1;

     e1_trials ( :, :, i) = e1 .* conj(e1);

     

     w2_trials ( :, :, i) = w2;

     e2_trials ( :, :, i) = e2 .* conj(e2);

end

w1_trials_ave = sum(w1_trials , 3) ./ trials;

e1_trials_ave = sum(e1_trials , 3) ./ trials;

w2_trials_ave = sum(w2_trials , 3) ./ trials;

e2_trials_ave = sum(e2_trials , 3) ./ trials;



sJmin = sum(Jmin) / trials;                         % 维纳误差 平均

Jex1 = e1_trials_ave - sJmin;                       % 剩余均方误差 

Jex2 = e2_trials_ave - sJmin;                       % 剩余均方误差 

sum_eig_trials = sum(sum_eig) / trials;             % 平均特征值之和

Jexfin1 = mu1 * sJmin * (sum_eig_trials / (2 - mu1 * sum_eig_trials));

Jexfin2 = mu2 * sJmin * (sum_eig_trials / (2 - mu2 * sum_eig_trials));

M1 = Jexfin1 / sJmin;                               % 失调参数

M2 = Jexfin2 / sJmin;                               % 失调参数

disp(M1);

disp(M2);

%% 画滤波器的权向量和学习曲线

figure(1);                                  % 画权值迭代过程

plot(w1(1, :), 'b');                        % 单次估计 w1

hold on

plot(w1(2, :), 'r');

hold on

plot(w1_trials_ave(1, :), 'b');             % 估计平均 w1

hold on

plot(w1_trials_ave(2, :), 'r');

xlabel('迭代次数n');

ylabel('抽头权值');

title('滤波器权系数的更新过程 \mu=0.05');



figure(2);                                  % 画权值迭代过程

plot(w2(1, :), 'b');                        % 单次估计 w2

hold on

plot(w2(2, :), 'r');

hold on

plot(w2_trials_ave(1, :), 'b');             % 估计平均 w2

hold on

plot(w2_trials_ave(2, :), 'r');

xlabel('迭代次数n');

ylabel('抽头权值');

title('滤波器权系数的更新过程 \mu=0.005');



figure(3);                                  % 画学习曲线

plot(Jex1, 'b');                            % 画mu1时的学习曲线

hold on

plot(Jex2, 'r');                            % 画mu2时的学习曲线

hold on

axis([0,550,0,1]);

xlabel('迭代次数n');

ylabel('MSE');

title('LMS算法的学习曲线');

legend('\mu1=0.05', '\mu2=0.005');

#include"math.h"



#define PI 3.1415926

#define COEFFNUMBER 16

#define INPUTNUMBER 1024



int FIRLMS(int *nx,float *nh,int nError,int nCoeffNumber);



float h[COEFFNUMBER],fU;

int xx[INPUTNUMBER],rr[INPUTNUMBER],wc[INPUTNUMBER];



main()

{

        int i,nLastOutput;

        

        nLastOutput=0;

        fU=0.0005;

        for ( i=0;i<COEFFNUMBER;i++ )        h[i]=0;

        for ( i=0;i<INPUTNUMBER;i++ )

        {

                xx[i]=256*sin(i*2*PI/34);

                rr[i]=wc[i]=0;

        }

        for ( i=COEFFNUMBER+1;i<INPUTNUMBER;i++ )

        {

                nLastOutput=FIRLMS(xx+i,h,nLastOutput-xx[i-1],COEFFNUMBER);        // break point

                rr[i]=nLastOutput;

                wc[i]=rr[i]-xx[i];

        }

        exit(0);

}



int FIRLMS(int *nx,float *nh,int nError,int nCoeffNumber)

{

        int i,r;

        float fWork;

        

        r=0;

        for ( i=0;i<nCoeffNumber;i++ )

        {

                fWork=nx[i]*nError*fU;

                nh[i]+=fWork;

                r+=(nx[i-i]*nh[i]);

        }

        r/=128;

        return r;

}

// LMS FIR (delayed version)

ushort dlms (DATA *x, DATA *h, DATA *r, DATA *des,DATA *dbuffer, DATA step, ushort nh, ushort nx)

// Adaptive delayed LMS filter (fast implemented)

ushort oflag = dlmsfast (DATA *x, DATA *h, DATA *r, DATA*des, DATA *dbuffer, DATA step, ushort nh, ushort nx)

http://xiahouzuoxin.github.io/notes/html/Kalman%E6%BB%A4%E6%B3%A2%E5%99%A8%E4%BB%8E%E5%8E%9F%E7%90%86%E5%88%B0%E5%AE%9E%E7%8E%B0.html

clc; clear all; close all;

%% Necessary parameter

N = 4;                      % 状态向量维数

M = 1;                      % 观测向量维数

samples = 512;              % 样本量

trails = 100;               % 试验次数

den = [1, -1.6, 1.46, -0.616, 0.1525];   %AR parameter

num = 1;

noise = sqrt(0.0332) * randn(samples, trails);

signal = filter(num, den, noise);



%% Kalman Filter Algorithm

%initial

Q2 = 0.005;                 % 观察噪声自相关矩阵

signal_u = zeros(samples + 4, trails);     % 观测信号 u(n)

u = zeros(N, samples);                     % 

c = zeros(samples, N);                     % 观察矩阵

K = zeros(N, samples);                     % 卡尔曼增益

A = zeros(num, samples);                   % 新息过程自相关矩阵

alpha = zeros(num, samples);               % 新息过程向量

omega = zeros(N, samples);                 % 预测器权向量

omega_aver = zeros(N, samples, trails);    %

P = zeros(N, N, samples);                  % 状态估计误差自相关矩阵

P(:, :, 1) = eye(N);                       % 

F = eye(N);                                % 转移矩阵

J = zeros(samples, trails);                % 均方误差

for m = 1 : trails

    % 产生观测信号signal_u(:, m),注意这里加0主要是为了n=1时signal_u正确取数

    signal_u(:, m) = [0; 0; 0; 0; signal(:, m)];

 %   signal_u(:, m) = [zeros(1, N - 1), signal(:, m)];

    for n = 2 : samples

        u(:, n) = signal_u(n + 3 : -1 : n, m).';

        c(n, :) = u(:, n)';

        alpha(n) = conj(signal(n, m)) - c(n, :) * omega(:, n - 1);

        A(n) = c(n, :) * P(:, :, n - 1) * c(n, :)' + Q2;

        K(:, n) = P(:, :, n - 1) * c(n, :)' /A(n);

        omega(:, n) = omega(:, n - 1) + K(:, n) * alpha(n);

        P(:, :, n) = (eye(4) - K(:, n) * c(n, :)) * P(:, :, n - 1);

        J(n, m) = sum(diag(P(:, :, n)));

    end

    omega_aver(:, :, m) = omega; 

end

omega_ave = sum(omega_aver, 3)./trails;

J_ave = sum(J, 2)./trails;



%% Draw the picture of weight vector (1 times trial)

figure(1);

n = 1 : samples;

plot(n, omega(1, :),'LineWidth', 2);

hold on

plot(n, omega(2, :),':r','LineWidth', 2);

hold on

plot(n, omega(3, :),'--g','LineWidth', 2);

hold on

plot(n, omega(4, :), '-.k','LineWidth', 2);

title('Using Kalman Filter algorithm to estimate weighter vector (1 times trial)');

xlabel('iterative times n');

ylabel('weight vector');

legend('\omega1','\omega2','\omega3','\omega4',0);



%% Draw the picture of weight vector (100 times trial)

figure(2)

plot(n, omega_ave(1, :),'LineWidth', 2);

hold on

plot(n, omega_ave(2, :),':r','LineWidth', 2);

hold on

plot(n, omega_ave(3, :),'--g','LineWidth', 2);

hold on

plot(n, omega_ave(4, :), '-.k','LineWidth', 2);

title('Using Kalman Filter algorithm to estimate weighter vector (100 times trial)');

xlabel('iterative times n');

ylabel('weight vector');

legend('\omega1','\omega2','\omega3','\omega4',0);



%% Draw the picture of MSE(Mean Square Error) (100 times trial)

figure(3);                                       

plot(n, J_ave, 'LineWidth', 2); 

xlabel('iterative times n');

ylabel('MSE');

title('基于Kalman算法的滤波器权系数的学习曲线');

http://processors.wiki.ti.com/index.php/GSG:Debugging_projects#Displaying_images

#include <stdio.h>



#define IMAG_EWIDTH   80

#define IMAG_EHEIGHT  80



void CreatImage(unsigned char *pImage,int nWidth,int nHeight);

void ReadImage(unsigned char *pImage,char *cFileName,int nWidth,int nHeight);

void Reverse(unsigned char *pImage, unsigned char *dbTargetImage, int nWidth,int nHeight);



unsigned char dbImage[IMAG_EWIDTH * IMAG_EHEIGHT];

unsigned char dbTargetImage[IMAG_EWIDTH * IMAG_EHEIGHT];

        

int main(void)

{

        CreatImage(dbImage,  IMAG_EWIDTH, IMAG_EHEIGHT);

        ReadImage(dbImage, "..\\1.bmp", IMAG_EWIDTH, IMAG_EHEIGHT);

        Reverse(dbImage, dbTargetImage, IMAG_EWIDTH, IMAG_EHEIGHT);

        

        while(1);

}                



void CreatImage(unsigned char *pImage,int nWidth,int nHeight)

{

        int x, y;

        int tmp1, tmp2;

        unsigned char *pWork = pImage;



        tmp1 = 256/ 16;

        tmp2 = nHeight / 16;



        for (y=0; y < nHeight;y ++ )

        {

                for (x = 0;x < nWidth; x++, tmp1++ )

                {

                        (*pWork) = (y / tmp2) * tmp1;

                }

        }

}





void ReadImage(unsigned char *pImage, char *cFileName, int nWidth, int nHeight)

{

        int j;

        unsigned char *pWork;

        FILE *fp;



        if (fp = fopen(cFileName, "rb"))

        {

                fseek(fp, 1078L, SEEK_SET);

                pWork = pImage + (nHeight - 1) * nWidth;

                for (j = 0; j < nHeight; j++, pWork -= nWidth)

                {

                        fread(pWork,nWidth,1,fp);

                }

                fclose(fp);

        }

}



void Reverse(unsigned char *pImage, unsigned char *dbTargetImage, int nWidth,int nHeight)

{

        int i, j;

        unsigned char *pImg = dbImage;

        unsigned char *pImg1 = dbTargetImage;

        for (i = 0; i < nHeight; i++ )

        {

                for (j = 0; j < nWidth; j++, pImg++, pImg1++)

                {

                        (*pImg1) = (~(*pImg)) & 0x0ff;

                }

        }

}

#define IMAGEWIDTH   80  // 图片宽度

#define IMAGEHEIGHT  80  // 图片长度



extern unsigned char dbImage[IMAGEWIDTH * IMAGEHEIGHT];     // 原图

extern unsigned char dbTargetImage[IMAGEWIDTH * IMAGEHEIGHT]; // 目标图



int mi, mj, m_nWork1, m_nWork2;    // 中间变量

unsigned int m_nWork, *pWork;     // 中间变量

unsigned char *pImg1, *pImg2, *pImg3, *pImg;  // 图像指针

unsigned int x11, x12, x13, x21, x22, x23, x31, x32, x33; // 矩阵 3x3



void laplace(int nWidth, int nHeight)

{

        int i;

        pImg = dbTargetImage;      // 目标图第一行元素清零

        for (i = 0; i < IMAGEWIDTH; i++, pImg++)

        {

                (*pImg) = 0;  

        }

        (*pImg) = 0;

        pImg1 = dbImage;   // 因为矩阵为3x3 所以先取图像的前三行

        pImg2 = pImg1 + IMAGEWIDTH;

        pImg3 = pImg2 + IMAGEWIDTH;

        for (i = 2; i < nHeight; i++)

        {

                pImg++;

                x11 = (*pImg1); pImag1++; 

                x12 = (*pImg1); pImag1++;         

                x21 = (*pImg2); pImag2++; 

                x22 = (*pImg2); pImag2++; 

                x31 = (*pImg3); pImag3++; 

                x32 = (*pImg3); pImag3++; 



                for (mi = 2; mi < nWidth; mi++, pImg++, pImg1++, pImg2++, pImg3++)

                {

                        x13 = (*pImg1); 

                        x23 = (*pImg2); 

                        x33 = (*pImg3);

                        

                        m_nWork1 = x22 << 2;   // 等同于

                        m_nWork1 += x22;       // m_nWork1 =  x22 * 5

                        m_nWork2 = x12 + x21 + x23 + x32;



                        m_nWork1 -= m_nWork2;

                        if (m_nWork1 > 255)

                        {

                                m_nWork1 = 255;

                        }

                        else if (m_nWork1 < 0)

                        {

                                m_nWork1 = 0;

                        }

                        (*pImg) = m_nWork1;

                        x11 = x12;

                        x12 = x12;

                        x21 = x22;

                        x22 = x23;

                        x31 = x32;

                        x32 = x33;

                }

                (*pImg) = 0;

                pImg++;

        }

}

http://jiaoyu.3158.cn/20130105/n5595270595830.html

unsigned char IntToALaw(int nInput)  //alaw转换，语音编码

{

        int segment;

        unsigned int i, sign,quant;

        unsigned int absol, temp;

        int nOutput;

        unsigned char cOutput;

        absol = abs(nInput);

        temp = absol;

        sign = (nInput >= 0)?1:0;

        for ( i=0;i<12;i++ )

        {

                nOutput = temp&0x4000;

                if ( nOutput )        

                {

                        break;

                }

                temp<<=1;

        }

        if ( i >= 12 )        

        {

                nOutput = 0;

        }

        else if( i>=7 && i<12 )

        {

                quant = (absol>>4)&0x0F;

                segment = 0;

                segment <<= 4; 

                nOutput = segment+quant;

        }              

        else

        {

                segment = 7-i;        //segment代表4-6位的数值

                quant = (absol>>(segment+3))&0x0F;

                segment <<= 4;

                nOutput = segment+quant;

        }

        if(sign)

                nOutput ^= 0xD5;        //正数与之相异或得到编码后的字节

        else

                nOutput ^= 0x55;        //负数与之相异或得到编码后的字节

        cOutput=(unsigned char)nOutput;

        return cOutput;

}

int ALawToInt(unsigned char nInput)                //alaw反转换，语音解码

{

        int sign,segment;

        int temp,quant,nOutput;

        temp = nInput^0xD5;        //异或返回编码值

        sign = (temp&0x80)>>7;        //返回符号值

        segment = temp&0x70; 

        segment >>= 4;

        if(segment == 0)

        {

                segment = 1;

            quant = temp&0x0F;        //得到abcd位

                quant <<= segment;

                quant = quant|0x01;

        }

        else

        {

                quant = temp&0x0F;

                quant += 0x10;

                quant <<= 1;

                quant = quant|0x01;

                quant <<= segment-1;//移位返回

        }

        //quant <<= 3;

        if ( sign )

                nOutput = -quant;

        else

                nOutput = quant;

        return nOutput;

}



unsigned int G711ALawEncode(int nLeft, int nRight)  //八位转十六位数据,nLeft是左声道数据，nRight是右声道数据

{

        unsigned char cL, cR;

        unsigned int uWork;

        cL = IntToALaw(nLeft);

        cR = IntToALaw(nRight);

        uWork = cL; uWork <<= 8; uWork |= cR;        //左声道高位，右声道低位

        return(uWork);

}