白噪声(matlab)

%-----------------------------------------------------------------

% To test rand.m and to generate the white noise signal u(n) 

%                   with uniform distribution

% 产生均匀分布的随机白噪信号,并观察数据分布的直方图

%-----------------------------------------------------------------

clear;

N=50000;           

u=rand(1,N);       

u_mean=mean(u)     

power_u=var(u)     

subplot(211)       

plot(u(1:100));grid on;

ylabel('u(n)')     

xlabel('n')     

subplot(212)       

hist(u,50);grid on;         

ylabel('histogram of u(n)')    

思考：如何产生一均匀分布、均值为0，功率为0.01的白噪声信号u(n)

提示：u（n）~[a，b]上均匀分布，E（u）= (a+b)/2, 

% to generate the white noise signal u(n)  with uniform distribution

% and power p;

% 产生均匀分布的白噪信号，使均值为0，功率为p

%-----------------------------------------------------------------

clear;

p=0.01;

N=50000;

u=rand(1,N);

u=u-mean(u);

a=sqrt(12*p); 

u1=u*a;

power_u1=dot(u1,u1)/N  

subplot(211)

plot(u1(1:200));grid on;

ylabel('u(n)')     

xlabel('n')     

example2: 产生零均值功率0.1且服从高斯分布的白噪声

%-----------------------------------------------------------------

% to test randn.m and to generate the white noise signal u(n) 

%                   with Gaussian distribution and power p

% 产生高斯分布的白噪信号，使功率为p，并观察数据分布的直方图

%-----------------------------------------------------------------

clear;

p=0.1;

N=500000;

u=randn(1,N);

a=sqrt(p)

u=u*a; 

power_u=var(u)

subplot(211)

plot(u(1:200));grid on;

ylabel('u(n)');

xlabel('n')     

subplot(212)

hist(u,50);grid on; 

ylabel('histogram of u(n)');

example3: sinc信号

%-----------------------------------------------------------------

% to generate the sinc function.

% 产生一 sinc 函数；

%-----------------------------------------------------------------

clear;

n=200;

stept=4*pi/n;

t=-2*pi:stept:2*pi;

y=sinc(t);

plot(t,y,t,zeros(size(t)));

ylabel('sinc(t)');

xlabel('t=-2*pi~2*pi');

grid on;

example4: chirp信号

%-----------------------------------------------------------------

% to test chirp.m and to generate the chirp signal x(t) 

% 产生一 chirp 信号；

% chirp(T0,F0,T1,F1):

% T0: 信号的开始时间；    F0:信号在T0时的瞬时频率，单位为Hz；

% T1: 信号的结束时间；    F1:信号在T1时的瞬时频率，单位为Hz；

%-----------------------------------------------------------------

  clear;

  t=0:0.001:1;   

  x=chirp(t,0,1,125);

  plot(t,x);

  ylabel('x(t)')     

  xlabel('t')

%-------------------------------------------------------------------------

% to test specgram.m ：估计信号谱图（SFFT）

%specgram（x，Nfft， Fs，window，Noverlap）

%x-信号；Fs抽样频率（2），Nfft做FFT长度（256）, window窗函数（Hanning）

%Noverlap：估计功率谱时每一段叠合长度（0）

%-------------------------------------------------------------------------

clear;

t=0:0.001:1.024-.001; 

N=1024; 

% 得到两个Chirp 信号，并相加；

y1=chirp(t,0,1,350);

y2=chirp(t,350,1,0);

y=y1+y2;

subplot(211);plot(t,y1);

ylabel(' Chirp signal y1')

% 求两个Chirp 信号和的短时傅里叶变换；

[S,F,T]=specgram(y,127,1,hanning(127),126); 

subplot(212);

surf(T/1000,F,abs(S).^2)

view(-80,30);

shading flat;

colormap(cool);

xlabel('Time')

ylabel('Frequency')

zlabel('spectrogram')

还有diric信号（周期SINC）gauspuls(高斯信号)pulstran(脉冲串信号)tripuls三角波脉冲信号

example5:线性卷积

%-----------------------------------------------------------------

% to test conv.m 

% 计算两个序列的线性卷积；

%-----------------------------------------------------------------

clear;

N=5;

M=6;

L=N+M-1;

x=[1,2,3,4,5];

h=[6,2,3,6,4,2];

y=conv(x,h);

nx=0:N-1;

nh=0:M-1;

ny=0:L-1;

subplot(231);

stem(nx,x,'.k');xlabel('n');ylabel('x(n)');grid on;

subplot(232);

stem(nh,h,'.k');xlabel('n');ylabel('h(n)');grid on;

subplot(233);

stem(ny,y,'.k');xlabel('n');ylabel('y(n)');grid on;

思考：

设信号x(n)由正弦信号加均值为0白噪声所组成，正弦信号幅度为1，白噪声方差为1

SNR = 10LG(PS/PU) = -3dB,试分析信号；

%-----------------------------------------------------------------

% to test xcorr.m

% 求两个序列的互相关函数，或一个序列的自相关函数；

%-----------------------------------------------------------------

clear;

N=500;

p1=1;

p2=0.1;

f=1/8;

Mlag=50;

u=randn(1,N);

n=[0:N-1];

s=sin(2*pi*f*n);

% 混有高斯白噪的正弦信号的自相关

u1=u*sqrt(p1);

x1=u1(1:N)+s;

rx1=xcorr(x1,Mlag,'biased');

subplot(221);

plot(x1(1:Mlag));

xlabel('n');

ylabel('x1(n)');grid on;

subplot(223);

plot((-Mlag:Mlag),rx1);grid on;

xlabel('m');ylabel('rx1(m)');

% 高斯白噪功率由原来的p1减少为p2，再观察混合信号的自相关

u2=u*sqrt(p2);

x2=u2(1:N)+s;

rx2=xcorr(x2,Mlag,'biased');

subplot(222);

plot(x2(1:Mlag));

xlabel('n');ylabel('x2(n)');grid on;

subplot(224);

plot((-Mlag:Mlag),rx2);

grid on;xlabel('m');ylabel('rx2(m)');

用乘同余法产生(见光盘 FLch2bzsheg2.m)

1编程如下：

A=6;  x0=1;  M=255;  f=2;   N=100；  %初始化；

x0=1;  M=255;

for k=1: N           %乘同余法递推100次；

   x2=A*x0;        %分别用x2和x0表示xi+1和xi-1；

  x1=mod (x2,M);    %取x2存储器的数除以M的余数放x1（xi）中；

   v1=x1/256;       %将x1存储器中的数除以256得到小于1的随机数放v1中；

   v(:,k)=(v1-0.5 )*f;  %将v1中的数（）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化；     

   x0=x1;           % xi-1= xi；

   v0=v1;

end                 %递推100次结束；

v2=v                %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中；

k1=k;

%grapher             %以下是绘图程序；

k=1:k1;

plot(k,v,k,v,'r');

xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声')

② 程序运行结果如图2.6所示。

图2.6  采用MATLAB产生的(-1,+1)均匀分布的白噪声序列

③ 产生的(-1，1)均匀分布的白噪声序列

在程序运行结束后，产生的(-1，1)均匀分布的白噪声序列，直接从MATLAB的window界面中copy出来如下（v2中每行存6个随机数）：

v2 =

   -0.9531   -0.7188    0.6875   -0.8359   -0.0156    0.9219

    0.5703    0.4531   -0.2500   -0.4844    0.1016   -0.3672

    0.8047   -0.1328    0.2188    0.3359   -0.9531   -0.7188

    0.6875   -0.8359   -0.0156    0.9219    0.5703    0.4531

   -0.2500   -0.4844    0.1016   -0.3672    0.8047   -0.1328 

    0.2188    0.3359   -0.9531   -0.7188    0.6875   -0.8359

   -0.0156    0.9219    0.5703    0.4531   -0.2500   -0.4844

    0.1016   -0.3672    0.8047   -0.1328    0.2188    0.3359 

   -0.9531   -0.7188    0.6875   -0.8359   -0.0156    0.9219 

    0.5703    0.4531   -0.2500   -0.4844    0.1016   -0.3672  

    0.8047   -0.1328    0.2188    0.3359   -0.9531   -0.7188

    0.6875   -0.8359   -0.0156    0.9219    0.5703    0.4531

   -0.2500   -0.4844    0.1016   -0.3672    0.8047   -0.1328

    0.2188    0.3359   -0.9531   -0.7188    0.6875   -0.8359

   -0.0156    0.9219    0.5703    0.4531   -0.2500   -0.4844

    0.1016   -0.3672    0.8047   -0.1328    0.2188    0.3359

   -0.9531   -0.7188    0.6875   -0.8359

显然，只要在例2.2程序的初始化部分中给N=300，f=6，运行程序就可以得到如图2.3所示的（-3，3）的白噪声过程.

①编程如下：

 A=6;  x0=1;  M=255;  f=6;   N=300；  %初始化；

x0=1;  M=255;

for k=1: N           %乘同余法递推100次；

   x2=A*x0;        %分别用x2和x0表示xi+1和xi-1；

  x1=mod (x2,M);    %取x2存储器的数除以M的余数放x1（xi）中；

   v1=x1/256;       %将x1存储器中的数除以256得到小于1的随机数放v1中；

   v(:,k)=(v1-0.5 )*f;  %将v1中的数（）减去0.5再乘以存储器f中的系数，存放在矩阵存储器v的第k列中，v(:,k)表示行不变、列随递推循环次数变化；     

   x0=x1;           % xi-1= xi；

   v0=v1;

end                 %递推100次结束；

v2=v                %该语句后无‘；’，实现矩阵存储器v中随机数放在v2中，且可直接显示在MATLAB的window中；

k1=k;

%grapher             %以下是绘图程序；

k=1:k1;

plot(k,v,k,v,'r');

xlabel('k'), ylabel('v');tktle(' (-1,+1)均匀分布的白噪声')

② 程序运行结果如图2.3所示。

图2.3 白噪声过程