matlab_seven.m

随机信号的一些基本MATLAB程序

%说明：此程序文件为随机信号分析与处理中第七讲中程序，包括随机信号模拟、数字特征估计、概率密度估计、功率谱估计及其性能分析，共12个小程序
%编写：张文明
%程序1、产生均匀分布白噪声随机数矩阵
x=rand(1000,1);
subplot(121);
plot(x)
subplot(122);
x1=0:0.05:1;
hist(x,x1)

% 程序2、产生标准正态分布白噪声随机数矩阵
% x=randn(1000,1);
% subplot(121);
% plot(x)
% subplot(122);
% x1=-4:0.1:4;
% hist(x,x1)

%程序3、产生韦伯分布随机数矩阵
% x=weibrnd(1,1.5,300,1);
% subplot(121);
% plot(x)
% subplot(122);
% x1=0:0.1:4;
% hist(x,x1)

%程序4、相关正态随机序列的产生
clear;
a=0.8;
a1=0.01;
sigma=2;
N=500;
u=randn(N,1);
x(1)=sigma*u(1)/sqrt(1-a^2);
for i=2:N
    x(i)=a*x(i-1)+sigma*u(i);
end
subplot(121);
plot(x);
subplot(122);
y=xcorr(x);
y=y./max(y);
for i=1:41
    y1(i)=y(i+480);
end
k=-20:20;
y2=4/(1-0.8^2)*0.8.^abs(k);
y2=y2./max(y2);
plot(k,y1,k,y2);


%程序5、随机序列均值、均方值、均方根值、标准差、方差计算程序
% clear;
% N=10000;
% y=randn(1,N);%产生N个标准正态分布白噪声序列
% plot(y);
% yMean=mean(y);   %均值
% disp('均值');yMean   
% y2p=y*y'/N;      %均方值计算
% disp('均方值');y2p
% ysp=sqrt(y2p);   %均方根值计算
% disp('均方根值');ysp
% ystd=std(y,1);     %均方差（标准差）值计算
% disp('均方差（标准差）值');ystd
% yd=ystd.*ystd;    %方差计算
% disp('方差');yd


%程序6、随机序列自相关函数分析
% clf;
% N=1000;
% n=0:N-1;
% Fs=500;
% t=n/Fs;
% Lag=100;
% x=sin(2*pi*10*t)+0.6*randn(1,length(t));  %信号产生
% [c,lags]=xcorr(x,Lag,'unbiased');
% subplot(221)
% plot(t,x)
% xlabel('t');ylabel('x(t)');
% title('Original Signal x');
% grid
% subplot(222)
% plot(lags/Fs,c);
% xlabel('t');ylabel('Rx(t)');
% title('Autocorrelation');
% grid
% x1=randn(1,length(x));%信号产生
% [c,lags]=xcorr(x1,Lag,'unbiased');
% subplot(223)
% plot(t,x1)
% xlabel('t');ylabel('x1(t)');
% title('Original Signal x1');
% grid
% subplot(224)
% plot(lags/Fs,c);
% xlabel('t');ylabel('Rx1(t)');
% title('Autocorrelation');
% grid


%程序7、利用互相关函数计算时移程序
% clf
% N=1000;
% n=0:N-1;
% Fs=500;
% t=n/Fs;
% Lag=200;
% x1=90*sinc(pi*(n-0.1*Fs));%signal x1 y1
% y1=50*sinc(pi*(n-0.3*Fs));
% [c,lags]=xcorr(x1,y1,Lag,'unbiased');
% subplot(121)
% plot(t,x1,'r')
% hold on
% plot(t,y1,'b');
% xlabel('t');ylabel('x(1t) y1(t)');
% title('Original Signal ');
% grid
% hold off
% subplot(122)
% plot(lags/Fs,c,'r')
% xlabel('t');ylabel('Rxy(t)');
% title('Correlation ');
% grid

%程序8、估计概率密度程序
% a=0.8;
% sigma=2;
% N=200;
% u=randn(N,1);
% x(1)=sigma*u(1)/sqrt(1-a^2);
% for i=2:N
%     x(i)=a*x(i-1)+sigma*u(i);
% end
% [f,xi] = ksdensity(x); 
% subplot(121);
% plot(xi,f); 
% xlabel('x');
% ylabel('f(x)');
% axis([-15 15 0 0.4]);
% subplot(122);
% k=-15:15;
% hist(x,k); 
% xlabel('x');
% ylabel('f(x)');


%程序9、利用周期图法估计功率谱密度程序
% clf
% Fs=1000;
% N=256;   %case 1 :N=256
% Nfft=256;
% n=0:N-1;
% t=n/Fs;
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);
% f=(0:length(Pxx)-1)*Fs/length(Pxx);
% subplot(211)
% plot(f,Pxx);
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Periodogram N=256 ');
% grid
% 
% N=1024;   %case 2 :N=1024
% Nfft=1024;
% n=0:N-1;
% t=n/Fs;
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);
% f=(0:length(Pxx)-1)*Fs/length(Pxx);
% subplot(212)
% plot(f,Pxx);
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Periodogram N=1024 ');
% grid




%程序10、利用分段周期图法估计功率谱密度程序
% clf
% Fs=1000;
% N=1024;%case 1 :no overlap
% Nsec=256;
% n=0:N-1;
% t=n/Fs;
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;
% Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;
% Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;
% Pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;
% Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
% f=(0:length(Pxx)-1)*Fs/length(Pxx);
% subplot(211)
% plot(f,Pxx);
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Averaged Periodogram(no overlap) N=4*256 ');
% grid
% 
% N=1024;%case 2 :half overlap
% Nsec=256;
% n=0:N-1;
% t=n/Fs;
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec;
% Pxx2=abs(fft(xn(129:384),Nsec).^2)/Nsec;
% Pxx3=abs(fft(xn(257:512),Nsec).^2)/Nsec;
% Pxx4=abs(fft(xn(385:640),Nsec).^2)/Nsec;
% Pxx5=abs(fft(xn(513:768),Nsec).^2)/Nsec;
% Pxx6=abs(fft(xn(641:896),Nsec).^2)/Nsec;
% Pxx7=abs(fft(xn(769:1024),Nsec).^2)/Nsec;
% Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
% f=(0:length(Pxx)-1)*Fs/length(Pxx);
% subplot(212)
% plot(f,Pxx);
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Averaged Periodogram(half overlap) N=1024 ');
% grid


%程序11、利用加窗周期图法估计功率谱密度程序
% clf
% Fs=1000;
% N=1024;%case 1 :no overlap
% Nsec=256;
% n=0:N-1;
% t=n/Fs;
% w=hanning(256)';
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% Pxx1=abs(fft(w.*xn(1:256),Nsec).^2)/Nsec;
% Pxx2=abs(fft(w.*xn(257:512),Nsec).^2)/Nsec;
% Pxx3=abs(fft(w.*xn(513:768),Nsec).^2)/Nsec;
% Pxx4=abs(fft(w.*xn(769:1024),Nsec).^2)/Nsec;
% Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);
% f=(0:length(Pxx)-1)*Fs/length(Pxx);
% subplot(211)
% plot(f,Pxx);
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Averaged Modified Periodogram(no overlap) N=4*256 ');
% grid
% 
% N=1024;%case 2 :half overlap
% Nsec=256;
% n=0:N-1;
% t=n/Fs;
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% Pxx1=abs(fft(w.*xn(1:256),Nsec).^2)/Nsec;
% Pxx2=abs(fft(w.*xn(129:384),Nsec).^2)/Nsec;
% Pxx3=abs(fft(w.*xn(257:512),Nsec).^2)/Nsec;
% Pxx4=abs(fft(w.*xn(385:640),Nsec).^2)/Nsec;
% Pxx5=abs(fft(w.*xn(513:768),Nsec).^2)/Nsec;
% Pxx6=abs(fft(w.*xn(641:896),Nsec).^2)/Nsec;
% Pxx7=abs(fft(w.*xn(769:1024),Nsec).^2)/Nsec;
% Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
% f=(0:length(Pxx)-1)*Fs/length(Pxx);
% subplot(212)
% plot(f,Pxx);
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Averaged Modified Periodogram(half overlap) N=1024 ');
% grid

%程序12、利用最大熵法与加窗改进周期图法估计功率谱密度程序
% clf
% Fs=1000;
% N=1024;%case 1 :no overlap
% Nfft=256;
% n=0:N-1;
% t=n/Fs;
% randn('state',0);
% xn=sin(2*pi*50*t)+2*sin(2*pi*60*t)+randn(1,N);%signal xn
% [Pxx,f]=pmem(xn,54,Nfft,Fs);
% subplot(211)
% plot(f,10*log10(Pxx));
% xlabel('Frequency (Hz)');ylabel('Power Spectrum (dB)');
% title('Maxmum Entropy Method(MEM) Order=74');
% grid
% 
% window=hanning(256);
% noverlap=128;
% dflag='none';
% subplot(212)
% psd(xn,Nfft,Fs,window,noverlap,dflag);
% grid on