随机信号分析大作业——平稳随机过程的抽取插值方法的探讨

京嘟嘟

已于 2023-02-01 17:26:54 修改

阅读量2.5k

点赞数 3

分类专栏： MATLAB数字信号处理基础文章标签： matlab 开发语言 Powered by 金山文档

于 2023-01-31 17:58:51 首次发布

本文链接：https://blog.csdn.net/Indigo1234/article/details/128819715

版权

MATLAB数字信号处理基础专栏收录该内容

1 篇文章 0 订阅

订阅专栏

本实验代码参照了网上的相关代码，并进行了大量的修改和补充。注释版代码我会放在文后。

一、实验目的

（1）了解确定信号的采样与平稳随机信号的采样之间的关系，掌握信号的采样定理及其应用；

（2）掌握随机信号的均值、方差、自相关函数、概率密度、频谱及功率谱密度的特性；

（3）掌握随机信号的分析方法；

（4）熟悉常用的信号处理仿真软件平台：MATLAB或C/C++。

二、实验内容

（一）实验原理

确定信号的采样符合香农定理，那么随机信号的采样是否符合香农定理呢？答案是肯定的。香农定理可以推广到随机信号的采样。

若X(t)为平稳随机过程，且具有零均值，它的功率谱密度 $\text{[math]}$ 限于 $\text{[math]}$ 之间。当满足条件 $\text{[math]}$ 时，便可将 $\text{[math]}$ 按照它的振幅样本展开为

$\text{[math]}$

上式为平稳随机过程的采样定理，式中T为采样周期。

（二）实验任务与要求

①通过实验掌握随机信号的抽取与插值的基本原理、方法及其在实际应用中的意义。计算信号的数字特征，理解它们的物理概念。用MATLAB或C/C++软件仿真之一编写程序和仿真。系统框图如图6-2-1所示。

图6-2-1 抽取、插值系统框图

②输入信号 $\text{[math]}$ ：是在一个波形周期有80个采样点的正弦信号 $\text{[math]}$ 。信号的真实频率取决于采样频率。假设采样频率为8000Hz，则输入信号频率就是100Hz。 $\text{[math]}$ 为高斯白噪声。输入信号 $\text{[math]}$ 的时域、频域波形图如图6-2-2所示，高斯白噪声的时域、频域波形图如图6-2-3所示。

要求测试白噪声的的均值、均方值、方差、相关函数、概率密度、频谱及功率谱密度，并用波形图表示。分析实验结果，掌握均值、均方值、方差、自相关函数、频谱及功率谱密度的物理意义。

③滤波器设计：要求信号经滤波器后保留有用信号。滤波器的结构及参数和所涉及的采样频率取值根据滤波器在这里所起的作用、输入信号本身的特点所确定，在这里不作规定。

滤波器设计好之后要求测试它的频率特性，并将频率特性曲线用波形图来表示，以验证是否符合要求。

④信号抽取：所谓信号抽取也就是信号采样率的降低。将低通滤波器输出信号作为原始信号，用抽取的方法分别为以M = 2、M = 4对原始信号进行抽取。

⑤信号插值：将M = 2、M = 4抽取信号分别作为原始信号，用插值的方法，分别还原原始信号（滤波后的信号）。对原始信号进行抽取

⑥计算 $\text{[math]}$ 、 $\text{[math]}$ 、 $\text{[math]}$ 、 $\text{[math]}$ 的均值、均方值、方差、自相关函数、频谱、频率及功率谱密度。

⑦对采样前后、插值前插值后信号进行比较。观察信号频谱的变化及周期延拓性。

⑧讨论x(t)的自相关函数、功率谱密度与y(t)的自相关函数、功率谱密度之间的关系。

三、代码部分

%% 生成输入信号x(t)及相应的加性高斯白噪声
Fs = 40000;                                            % 生成模拟信号采样点的频率设为40000Hz
fs = 1/Fs;
Fx = 100;                                               % 输入信号频率
fx = 1/Fx; 
n = 5000;                                               % 设总共采样5000个点，一个周期80个采样点

x = 0: n-1;                                            % 得到离散序列
t = x/Fs;                                               % 得到时间采样序列
f = Fs*x/n ;                                           % 频域横坐标
nt = randn(1, n);                                   % 生成高斯随机信号
xt = sin(t*Fx*2*pi);                              % 生成原始信号
Xt.sig = nt + xt;                                    % 生成实际信号（原始信号带噪声）

% 绘制图像
figure(1);
subplot(3, 1, 1);
plot(xt); title("原始信号波形"); xlabel("t");
subplot(3, 1, 2);
plot(nt); title("产生的高斯信号");
subplot(3, 1, 3);
plot(Xt.sig); title("输入的实际信号");

%% 测试白噪声的均值、均方值、方差等参数
% 白噪声（即n(t)）
Nt.sig = nt;
Nt.mean = mean(Nt.sig);                                     % 求均值
Nt.var = var(Nt.sig, 1);                                    % 求方差
Nt.quad = Nt.var + Nt.mean.^2;                           % 求均方值，其等于方差加上均值的平方
Nt.acor = xcorr(Nt.sig);                                   % 求自相关函数
Nt.ksde = ksdensity(Nt.sig);                             % 求概率密度
Nt.fft = fft(Nt.sig);                                        % 求频谱
% Xt.psd = abs(Xt.fft).^2/(length(Xt.fft)*fs);  
Nt.psd = Nt.fft.*conj(Nt.fft)/n;                       % 求功率谱密度 

% 绘制图像
figure(31);
subplot(2, 2, 1);
plot(Nt.acor); 
title("n(t)的自相关函数");

subplot(2, 2, 2);
ksdensity(Nt.sig); 
title("n(t)的概率密度");

subplot(2, 2, 3);
plot(f(1: 300), abs(Nt.fft(1: 300))); 
title("n(t)的频谱");

subplot(2, 2, 4);
plot(Nt.psd(1:500));
title("n(t)的功率谱密度");

% 合成后的加性高斯白噪声（即Xt）
Xt.mean = mean(Xt.sig);                                     % 求均值
Xt.var = var(Xt.sig, 1);                                    % 求方差
Xt.quad = Xt.var + Xt.mean.^2;                           % 求均方值，其等于方差加上均值的平方
Xt.acor = xcorr(Xt.sig);                                   % 求自相关函数
Xt.ksde = ksdensity(Xt.sig);                             % 求概率密度
Xt.fft = fft(Xt.sig);                                        % 求频谱
% Xt.psd = abs(Xt.fft).^2/(length(Xt.fft)*fs);  
Xt.psd = Xt.fft.*conj(Xt.fft)/n;                       % 求功率谱密度 

% 绘制图像
figure(30);
subplot(2, 2, 1);
plot(Xt.acor); 
title("X(t)的自相关函数");

subplot(2, 2, 2);
ksdensity(Xt.sig); 
title("X(t)的概率密度");

subplot(2, 2, 3);
plot(f(1: 300), abs(Xt.fft(1: 300))); 
title("X(t)的频谱");

subplot(2, 2, 4);
plot(Xt.psd(1:500));
title("X(t)的功率谱密度");

%% 滤波器设计一（使用低通数字滤波器）
% 滤波器设置参数为带通，保留原始输入正弦信号（频率100Hz）
fp100 = 1000; fs100 = 2000;                                 % 载波频率为250Hz的条幅信号采取低通滤波器，其通带边界（截止）频率和阻带截止频率        
alphap100 = 1; alphas100 = 60;    

wp100 = 2*pi*fp100/Fs; ws100 = 2*pi*fs100/Fs; 
wap = tan(wp100/2); was = tan(ws100/2); Fs1=Fs/Fs;

[N,Wn]= buttord(wap, was, alphap100, alphas100,'s') ;[z, p, k]= buttap(N);
[bp, ap]= zp2tf(z, p, k);%得到传输函数
[bs, as]= lp2lp(bp, ap,wap);%低通到低通，频谱变换
[bz, az]= bilinear(bs, as,Fs1/2);%将模拟滤波器传递函数转换为数字滤波器传递函数
[H, w]= freqz(bz, az,256,Fs1*10000);%特性分析

% 绘制图像
figure(2);
plot (w, abs (H));% 画出数字滤波器的波形图
title('低通滤波器的频谱特性');
xlabel('f/hz');
grid on;

%% 将信号进行滤波，并求均值、均方值、方差等参数
Yt.sig = filter(bz, az, Xt.sig);
Yt.mean = mean(Yt.sig);                                     % 求均值
Yt.var = var(Yt.sig, 1);                                    % 求方差
Yt.quad = Yt.var + Yt.mean.^2;                           % 求均方值，均方值表示信号的平均功率
Yt.acor = xcorr(Yt.sig);                                   % 求自相关
Yt.ksde = ksdensity(Yt.sig);                             % 求概率密度
Yt.fft = fft(Yt.sig);                                        % 求频谱
% Yt.psd = abs(Yt.fft).^2/(length(Yt.fft)*fs);     % 求功率谱密度 
Yt.psd = Yt.fft.*conj(Yt.fft)/n;

% 绘制图像
figure(3);
subplot(2, 1, 1);
plot(Xt.sig); 
title("x(t)滤波前图像");
subplot(2, 1, 2);
plot(Yt.sig); 
title("y(t)滤波后图像");

figure(4);
subplot(2, 2, 1);
plot(Yt.acor); 
title("y(t)的自相关函数");

subplot(2, 2, 2);
ksdensity(Yt.sig); 
title("y(t)的概率密度");

subplot(2, 2, 3);
plot(f(1: 300), abs(Yt.fft(1: 300))); 
title("y(t)的频谱");

subplot(2, 2, 4);
plot(Yt.psd(1:500));
title("y(t)的功率谱密度");

%% 对实际信号处理，得到去除相邻采样点间零点处的采样信号，为计算频谱和功率谱密度准备
Fs1=8000;                                                          % 对输入的实际信号采样周期设为8000Hz
t2=(0: Fs/Fs1: n-1)/Fs;                                      % 得到采样信号横坐标
sample=zeros(1,length(x)) ;                               % 得到采样信号横坐标向量
for k = 1 : length(x)
    if mod (k, (Fs/Fs1)) == 1
        sample(k) = Yt.sig(k);
    else
        sample(k) = 0;
    end
end

sample1 = zeros (1, ceil(length(x)/(Fs/Fs1)) ) ;
i = 1;
for k = 1: length(x)
    if mod(k, (Fs/Fs1)) == 1
        sample1(i) = Yt.sig(k); i = i+1;
    end
end

% 绘制图像
figure(5) ;
stem(sample1);% 除去相邻采样点间零点的采样信号
xlabel('t'); ylabel('y');
title('对实际信号采样后的信号');

%% 采样后信号xt（采样率为8000Hz）
Sxt.sig = sample1;                                                % 得到采样后的信号xt
Sxt.mean = mean(Sxt.sig);                                     % 求均值
Sxt.var = var(Sxt.sig, 1);                                    % 求方差
Sxt.quad = Sxt.var + Sxt.mean.^2;                          % 求均方值，均方值表示信号的平均功率
Sxt.fft = fft(sample1);                                        % 求频谱
Sxt.acor = xcorr(sample1);                                    % 求自相关
Sxt.psd = Sxt.fft.*conj(Sxt.fft)/n ;                     % 求功率谱密度PSD

% t1=(-n+1:n-1)/Fs ;                                             % 相关函数图像横坐标

% 绘制图像
figure (6);
subplot(3, 1, 1);
plot(abs(Sxt.fft(1 :800)));
xlabel('f') , ylabel('Y');
title('采样后信号频谱');

subplot(3, 1, 2);
plot(Sxt.acor);
xlabel('t'), ylabel('R');
title('采样后信号相关函数');

subplot(3, 1, 3);
plot(abs(Sxt.psd(1 : 500)));
xlabel('f'), ylabel('S');
title('采样后信号功率谱密度');

%% 间隔两点抽样（m = 2)
sample4 = zeros (1, ceil(length(t)/(Fs/Fs1))) ;     % 分配存储空间

for k=1:length(t2)                                               % 得到间隔两点的抽样信号
    if mod (k,3)==1
        sample4(k)=sample(k);
    else
        sample4(k)=0;
    end
end

sample5=zeros (1, ceil(length(t)/(Fs/Fs1)/3));     % 得到间隔两点的抽样信号（经过去零处理）
i=1;
for k=1:3:length(t2)
    sample5(i) = sample1(k);
    i=i+1;
end

% 绘制图像
figure (7);
plot(sample5);                                                      % 除去相邻采样点间零点的采样信号
xlabel('t'), ylabel('y');
title('（m = 2）抽取信号');

Sxt1.sig = sample4;                                               % 得到信号x1(t)
Sxt1.mean = mean(sample4);                                       % 求均值
Sxt1.var=var(sample4);                                          % 求方差
Sxt1.quad = Sxt1.var + Sxt1.mean.^2;                       % 求均方值，均方值表示信号的平均功率
Sxt1.acor = xcorr(sample4);                                    % 求自相关
Sxt1.fft=fft(sample4);                                          % 求频谱
Sxt1.psd = Sxt1.fft.*conj(Sxt1.fft)/n;                   % 求功率谱密度

% t22=(-n+1)/Fs*Fs1:(n-1)/Fs*Fs1;%相关函数图像横坐标

% 绘制图像
figure(8);
plot(Sxt1.acor);
xlabel('t'), ylabel('R');
title('（m = 2）抽取信号相关函数');

[F1,y1] = ksdensity(sample4);                                  % 求概率密度，y1为取值，F1为概率
figure(9);
plot(y1,F1);
xlabel('y'),ylabel('F');
title('（m = 2）抽取信号概率密度');

figure(10);
plot(abs(Sxt1.fft(1:800)));
xlabel('f'),ylabel('Y');
title('（m = 2）抽取信号频谱');

figure(11);
plot(abs(Sxt1.psd(1:800)));
xlabel('f'),ylabel('S');
title('（m = 2）抽取信号功率谱密度');

%% 求2插值（m = 2，x1(t)）

t3 = zeros(1,(length(t2)-1)*2+length(t2));              % 插值信号横坐标向量

for k=1:length(t3)                                                  % 构建插值向量
    if mod(k,3)==1
        t3(k)=(1/Fs1)*fix(k/3);
    end
    if mod(k,3)==2
        t3(k)=(1/Fs1)*fix(k/3)+(1/Fs1)/3;
    end
    if mod(k,3)==0
        t3(k)=(1/Fs1)*(k/3-1)+(2/Fs1)/3;
    end
end

% 构建插值信号（m = 2），并求相关参数
Ixt1.sig = interp1(t2,sample1,t3,'cubic');              % 线性三次内插，将sample1插入t2，返回t3（即m = 2的插值信号）
Ixt1.mean = mean(Ixt1.sig);                                     % 求均值
Ixt1.var = var(Ixt1.sig);                                        % 求方差
Ixt1.quad = Ixt1.var + Ixt1.mean.^2;                        % 求均方值
Ixt1.acor = xcorr(Ixt1.sig);                                   % 求自相关
Ixt1.fft = fft(Ixt1.sig);                                        % 求频谱
Ixt1.psd = Ixt1.fft.*conj(Ixt1.fft)/length(t2);      % 求功率谱密度

% 绘制图像
figure(12);
plot(Ixt1.sig);
xlabel('t'),ylabel('y');
title('（m = 2）插值信号');

figure(13);
plot(Ixt1.acor);
xlabel('t'),ylabel('R');
title('（m = 2）插值信号相关函数');

[F9,y9]=ksdensity(Ixt1.sig);                                   % 求概率密度，y9为取值，F9为概率
figure(14);
plot(y9,F9);
xlabel('y'),ylabel('F');
title('（m = 2）插值信号概率密度');

figure(15);
plot(abs(Ixt1.fft(1:300)));
xlabel('f'),ylabel('Y');
title('（m = 2）插值信号频谱');

figure(16);
plot(Ixt1.psd(1:300));
xlabel('f'),ylabel('S');
title('（m = 2）插值信号功率谱密度');

%% 间隔四点抽样（m = 4)
sample2 = zeros (1, ceil(length(t)/(Fs/Fs1))) ;

for k=1:length(t2)
    if mod (k,5)==1
        sample2(k)=sample(k);
    else
        sample2(k)=0;
    end
end

sample3=zeros (1, ceil(length(t) /(Fs/Fs1)/5));                 
i=1;
for k=1:5:length(t2)
    sample3(i)=sample1(k);
    i=i+1;
end

% 绘制图像
figure (17);
plot(sample3);                                                          % 除去相邻采样点间零点的采样信号
xlabel('t'), ylabel('y');
title('（m = 4）抽取信号');

Sxt2.sig = sample2;
Sxt2.mean=mean(sample2);                                           % 求均值
Sxt2.std=std(sample2);                                              % 求均方差
Sxt2.var=var(sample2);                                              % 求方差
Sxt2.acor = xcorr(sample2);                                        % 求自相关
Sxt2.fft=fft(sample2);                                              % 求频谱
Sxt2.quad=Sxt2.var+Sxt2.mean.^2;                               % 求均方值

% t3=(-n+1)/Fs*Fs1:(n-1)/Fs*Fs1;%相关函数图像横坐标

% 绘制图像
figure(18);
plot(Sxt2.acor);%相关函数
xlabel('t'), ylabel('R');
title('（m = 4）抽取信号相关函数');

[F2,y2]=ksdensity(sample2);                                     % 求概率密度，y9为取值，F9为概率
figure(19);
plot(y2,F2);%概率密度
xlabel('y'),ylabel('F');
title('（m = 4）抽取信号概率密度');

figure(20);
plot(f(1:800),abs(Sxt2.fft(1:800)));%频谱
xlabel('f'),ylabel('Y');
title('（m = 4）抽取信号频谱');

P2=Sxt2.fft.*conj(Sxt2.fft)/n;%频谱函数的平方求时间平均
figure(21);
plot(abs(P2(1:800)));%功率谱密度
xlabel('f'),ylabel('S');
title('（m = 4）抽取信号功率谱密度');

%% 求4插值（m = 4，x2(t)）
t4 = zeros(1,(length(t2)-1)*4+length(t2));%插值信号的横坐标向量
% 进行插值
for k=1:length(t4)
    if mod(k,5)==1
        t4(k)=(1/Fs1)*fix(k/5);
    end
    if mod(k,5)==2
        t4(k)=(1/Fs1)*fix(k/5)+(1/Fs1)/5;
    end
    if mod(k,5)==3
        t4(k)=(1/Fs1)*fix(k/5)+(2/Fs1)/5;
    end
    if mod(k,5)==4
        t4(k)=(1/Fs1)*fix(k/5)+(3/Fs1)/5;
    end
    if mod(k,5)==0
        t4(k)=(1/Fs1)*(k/5-1)+(4/Fs1)/5;
    end 
end

% 构建插值信号（m = 4），并求相关参数
Ixt2.sig = interp1(t2,sample1,t4,'cubic');
Ixt2.mean = mean(Ixt2.sig);
Ixt2.var = var(Ixt2.sig);
Ixt2.quad = Ixt2.var + Ixt2.mean.^2;
Ixt2.acor = xcorr(Ixt2.sig);
Ixt2.fft = fft(Ixt2.sig);
Ixt2.psd = Ixt2.fft.*conj(Ixt2.fft)/length(t4);

% 绘制图像
figure(22);
plot(t4, Ixt2.sig);
xlabel('t'),ylabel('y');
title('（m = 4）插值信号');

% t5=(-Ns+1:(Fs/fs1)/5:Ns-1)/Fs;%相关函数图像横坐标

figure(23);
plot(Ixt2.acor);
xlabel('t'),ylabel('R');
title('（m = 4）插值信号相关函数');

[F5,y5]=ksdensity(Ixt2.sig);                                    % 求概率密度，y5为取值，F5为概率
figure(24);
plot(y5,F5);
xlabel('y'),ylabel('F');
title('（m = 4）插值信号概率密度');

% f1=5*fs1*(0:length(t4)-1)/length(t4);%频域横坐标

figure(25);
plot(abs(Ixt2.fft(1:300)));
xlabel('f'),ylabel('Y');
title('（m = 4）插值信号频谱');

figure(26);
plot(Ixt2.psd(1:300));
xlabel('f'),ylabel('S');
title('（m = 4）插值信号功率谱密度');

%% figures

figure(27); % 采样前后
subplot(3, 1, 1);
plot(Sxt.sig(1: 200));
title('采样信号波形');

subplot(3, 1, 2);% m = 2
plot(Sxt1.sig);
title('（m = 2）抽取信号波形');

subplot(3, 1, 3);% m = 4
plot(Sxt2.sig);
title('（m = 4）抽取信号波形');

%%
figure(28); % 插值前后
subplot(2, 2, 1);
plot(abs(Sxt1.fft));
title('（m = 2）插值前信号频谱');

subplot(2, 2, 2);
plot(abs(Sxt2.fft));
title('（m = 4）插值前信号频谱');

subplot(2, 2, 3);
plot(abs(Ixt1.fft(1:100)));
title('（m = 2）插值后信号频谱');

subplot(2, 2, 4);
plot(abs(Ixt2.fft(1:100)));
title('（m = 4）插值后信号频谱');

%%
figure(29);% xt与yt自相关函数、功率谱密度
subplot(3, 2, 1);
plot(Xt.acor);
title('滤波前信号自相关函数');

subplot(3, 2, 2);
plot(Xt.psd(1: 50));
title('滤波前信号功率谱密度');

subplot(3, 2, 3);
plot(Yt.acor);
title('滤波后信号自相关函数');

subplot(3, 2, 4);
plot(Yt.psd(1: 50));
title('滤波后信号功率谱密度');

subplot(3, 2, 5);
plot(Xt.acor-Yt.acor);
title('滤波前后信号自相关函数作差');

subplot(3, 2, 6);
plot(Xt.psd-Yt.psd);
title('滤波前后信号功率谱密度作差');