基于Scilab的数字信号处理（1）

进行数字信号处理时，首先需要会从数据文件中读取数据，处理后还要能把数据保存到文件中。
涉及到数据文件读取的函数有write,read,save,load。

首先查看write的帮助文档。
write(file-desc,a,[format])
write(file-desc,a,k,format)
这里file-desc为字符串，给出保存的文件名，a是一个实矩阵或字符串列向量，k是一个整数向量，format指定一个'Fortran'风格的格式，

必须用括号括起来，且不能与浮点数，整数或字符编辑模式混合，没有学过fortran不知道怎么弄，从网上找了个表格，附在文后。

这里file-desc为字符串，给出保存的文件名，a是一个实矩阵或字符串列向量，k是一个整数向量，format指定一个'Fortran'风格的格式，

必须用括号括起来，且不能与浮点数，整数或字符编辑模式混合。%io(2)代表的文件是屏幕，而%io(1)代表的文件是键盘。

read帮助文档

[x]=read(file-desc,m,n,[format])
[x]=read(file-desc,m,n,k,format)

file-desc为要读取的数据文件名，m，n指定读取的矩阵大小，如果不知道行数，可以设置m=-1，则把整个文件的数据都读进来。k为一个

整数或一个整数向量。format与write的用法一致。

实例1:

-->A=rand(3,2)                                          建立一个随机矩阵
 A  =
 
    0.2320748    0.8833888 
    0.2312237    0.6525135 
    0.2164633    0.3076091 

 

-->write('test1',A,'(2E13.7)');                    把矩阵写入文件test1

-->A=read('test1',3,2,'(2E13.7)');             把数据读入变量A，需要指定行和列
 
-->A
 A  =
 
    0.2320748    0.8833888 
    0.2312237    0.6525135 
    0.2164633    0.3076091 

 

-->B=read('test2',-1,2,'(2E13.7)');            读数据时忽略行数

-->B
 B  =
 
    0.2320748    0.8833888 
    0.2312237    0.6525135 
    0.2164633    0.3076091 

 

save(filename [,x1,x2,...,xn])
save(fd [,x1,x2,...,xn])
load(filename [,x1,...,xn])
load(fd [,x1,...,xn])
save和load分别用于把变量列表中的变量保存到文件中或从文件中读取，如果省略其中的变量列表，则会保存工作空间中的所有变量或回

复之前保存的所有变量。

scilab提供了随机数发生器rand()，它提供了多个工作模式并可以根据需要设置或返回随机分布类型及种子。一种使用方法如下，产生一个

2x2的高斯分布矩阵：
实例2：
-->rand(2,2,'gauss')                         高斯分布即期望为0，方差为1的正态分布
 ans  =
 
  - 0.1890526    0.4249849 
    0.7088748  - 0.1481575 

 

此外，新版本的Scilab还提供了一个功能更强大的随机数发生器grand。它可以有更多的分布类型，并可以设置分布类型中的参数，而rand

只提供高斯分布和均匀分布。因此采用这个函数更为灵活。
实例3：
-->grand(2,2,'nor',1,2)                      产生一个正态分布，且期望为1，方差为2的2x2矩阵
 ans  =
 
    1.8777722    1.3958656 
  - 1.0576091    1.7667825 

为了得到特定均值和协方差的随机矢量，可以对高斯分布的结果进行线性变换。比如为了得到具有N（M，C）分布的随机矢量，M为均值向

量，C为协方差矩阵，为正定对称矩阵，因此可以表示为C=L^TL，L为上三角矩阵。可以对分布为N（0,1）的随机矢量x进行变换Y=C^1/2x+M=Lx+M。L可以通过对称正定矩阵的Cholesky 分解得到，在scilab中可以用函数chol函数计算得。

实例4：
-->m=[-2;1;10]                                定义均值矢量
 m  =
 
  - 2.  
    1.  
    10. 
 
-->L=[3 2 1;2 3 2;1 2 3]               定义协方差矩阵，为对称正定矩阵
 L  =
 
    3.    2.    1. 
    2.    3.    2. 
    1.    2.    3. 
 -->L2=chol(L)                                求解协方差矩阵的Cholesky分解上三角矩阵

 L2  =
 
    1.7320508    1.1547005    0.5773503 
    0.           1.2909944    1.0327956 
    0.           0.           1.2649111 
 -->x=grand(3,1,'nor',0,1)            产生高斯分布的随机矢量
 x  =
 
    0.1788246 
  - 0.7169367 
    1.0433    
 
-->y=L2'*x+m                              对高斯分布的随机矢量进行线性变换，得到N（M，C）分布的随机矢量
 y  =
 
  - 1.6902668 
    0.2809275 
    10.682477 

 

 poly用于定义多项式，在scilab中，多项式也是一个数据类型，它通过poly函数进行定义，并且可以对其进行运算。其形式为

p=poly(a,vname, ["flag"])

a可以为向量或矩阵，如果给出的是向量，则定义多项式，如果给出的是矩阵，则定义特征多项式。定义的形式由["flag"]确定，

flag可以为roots或coeff。而vname为多项式中的变量符号。roots函数用于求多项式的根。

实例5：

-->a=poly([1 2 3],'x','coeff')
 a  =
 
               2  
    1 + 2x + 3x   
 
-->b=poly([1 2 3],'x','roots')
 b  =
 
                2   3  
  - 6 + 11x - 6x + x   
 
-->c=poly([1 2;3 4],'x','roots')
 c  =
 
              2  
  - 2 - 5x + x   
 
-->a+b+c
 ans  =
 
               2   3  
  - 7 + 8x - 2x + x   

-->a*b
 ans  =
 
              2     3     4    5  
  - 6 - x - 2x + 22x - 16x + 3x   
 -->a/b
 ans  =
 
                 2     
      1 + 2x + 3x      
    ---------------    
                2   3  
  - 6 + 11x - 6x + x   

在scilab中，有理式采用list数据格式进行存取，它包含四个基本元素，第一个为'r'，表示这是一个有理式，第二第三个元素分别为有理式的分

子和分母，而最后一个元素为[]或字符串。

 freq函数用于计算多项式在特定点上的值，还经常用于计算滤波器或系统的频域相应，其实也就是给出传输函数的形式，然后计算在各个点

上的值。它有两个形式

   
   

    
    

     
     
[x]=freq(A,B,C [,D],f)
[x]=freq(NUM,DEN,f)

     
     
A,B,C,D，NUM,DEC均用于确定系统传输函数，而f列向量，用于指定要求那些点的值。

     
     
实例6：

     
     
-->a
 a  =
 
               2  
    1 + 2x + 3x   
 
-->freq(a,1,[1;2])
 ans  =
 
    6.    17. 

     
     
线性系统传输函数可以用有理式或状态空间法进行描述，在scilab中，采用syslin函数定义线性系统函数。
其格式为:
[sl]=syslin(dom,A,B,C [,D [,x0] ])
[sl]=syslin(dom,N,D)
[sl]=syslin(dom,H)
dom取值为‘c’，‘d’或[]或一个标量，表示系统为连续，离散或其它。A,B,C,D用于在状态空间中定义系统。
x0为矢量，给出了初始状态
N/D为多项式矩阵
H为有理矩阵或线性状态空间表示
而sl表示线性系统，很多函数可以接受这样的数据类型，其实它也是一个list。

     
     
实例6：

     
     
 a  =
 
               2  
    1 + 2x + 3x   
 
-->b
 b  =
 
                2   3  
  - 6 + 11x - 6x + x   
 
-->syslin('d',a,b)
 ans  =
 
                 2     
      1 + 2z + 3z      
    ---------------    
                2   3  
  - 6 + 11z - 6z + z   
 scilab还提供了系统并联，串联以及传输函数和状态空间表示之间的变换，并且还可以在连续时间系统和离散时间系统之间进行转换。

    
    

   
   

 

表一  可重复编辑符

 

表二 不可重复编辑符