精华内容
下载资源
问答
  • Matlab 产生白噪声和有色噪声序列

    万次阅读 2015-05-01 08:56:53
    一、白噪声和有色噪声定义 1.白噪声(white noise) 系统辨识中所用到的数据通常都是含有噪声的。从工程实际出发,这种噪声往往可以视为具有有理谱密度的平稳随机过程。白噪声是一种最简单的随机过程,是有一...

    原文地址 http://blog.sina.com.cn/s/blog_84024a4a01016fmb.html

    一、白噪声和有色噪声定义

    1.白噪声(white noise)

    系统辨识中所用到的数据通常都是含有噪声的。从工程实际出发,这种噪声往往可以视为具有有理谱密度的平稳随机过程。白噪声是一种最简单的随机过程,是有一系列不相关的随机变量组成的理想化随机过程。其自相关函数为dirac函数。

    2.有色噪声(colored noise)

    理想的白噪声只是一种理论上的抽象,在物理上是很难实现的,现实中并不存在这样的噪声。因而,工程实际中测量数据所包含的噪声往往是有色造势。所谓有色噪声(或相关噪声)是指序列中没一时刻的噪声相关。有色噪声可以看成是由白噪声序列驱动的线性环节的输出。

    二、白噪声与有色噪声区别

    (1)其实由定义可以看出,白噪声不同时刻是不相关的,自相关函数为脉冲函数;有色噪声则是相关的。

    (2)实际测试可以通过测试功率谱来区别,白噪声的功率谱在各频率的值都比较平均,有色噪声则会有较为明显的峰值。

    三、具体实例

    1.产生有色噪声e(k) = x(k) + 0.5*x(k-1)。其中,x(k)为方差为1的白噪声

    clear all; close all;
    clc
    L=500;  %仿真长度
    c = [1 -0.5];
    nc = length(c) - 1;
    xik=zeros(nc,1);  %白噪声初值
    xi=randn(L,1);  %产生均值为0,方差为1的高斯白噪声序列

    for k=1:L
        e(k)=c*[xi(k);xik];  %产生有色噪声
        %数据更新
        for i=nc:-1:2
            xik(i)=xik(i-1);
        end
        xik(1)=xi(k);
    end

    subplot(2,1,1);
    plot(xi);
    xlabel('k');ylabel('噪声幅值');title('白噪声序列');
    subplot(2,1,2);
    plot(e);
    xlabel('k');ylabel('噪声幅值');title('有色噪声序列');

    %测试功率谱

    [y1,f1] = Spectrum_Calc(xi',512);
    p1 = 1/L * y1.*conj(y1);

    figure(2)
    subplot(211)
    plot(f1,p1)

    [y2,f2] = Spectrum_Calc(e,512);
    p2 = 1/L * y2.*conj(y2);
    subplot(212)
    plot(f2,p2)

     

    运行结果:

    image

    image

    2:

     

    image

    clear all; close all;
    L=500;  %仿真长度
    d=[1 -1.5 0.7 0.1]; c=[1 0.5 0.2];  % 分子分母多项式系数
    nd=length(d)-1 ;nc=length(c)-1;   %阶次
    xik=zeros(nc,1);  %白噪声初值
    ek=zeros(nd,1);
    xi=randn(L,1);  %产生均值为0,方差为1的高斯白噪声序列

    for k=1:L
        e(k)=-d(2:nd+1)*ek+c*[xi(k);xik];  %产生有色噪声
        %数据更新
        for i=nd:-1:2
            ek(i)=ek(i-1);
        end
        ek(1)=e(k);
        for i=nc:-1:2
            xik(i)=xik(i-1);
        end
        xik(1)=xi(k);
    end
    subplot(2,1,1);
    plot(xi);
    xlabel('k');ylabel('噪声幅值');title('白噪声序列');
    subplot(2,1,2);
    plot(e);
    xlabel('k');ylabel('噪声幅值');title('有色噪声序列');

    %测试功率谱

    [y1,f1] = Spectrum_Calc(xi',512);
    p1 = 1/L * y1.*conj(y1);

    figure(2)
    subplot(211)
    plot(f1,p1)

    [y2,f2] = Spectrum_Calc(e,512);
    p2 = 1/L * y2.*conj(y2);
    subplot(212)
    plot(f2,p2)

    image

    image

    这个功率谱比例1效果要好很多。

    四、程序及结果说明:

    1.Spectrum_Calc函数为快速傅里叶变换函数,具体见博文

    http://blog.sina.com.cn/s/blog_84024a4a01015rez.html

    2.例2中提供了一个线性系统,将系统转换为差分方程形式可能更有助于程序的理解。系统的等价差分方程为:

    e(k) – 1.5e(k-1) + 0.7e(k-2) + 0.1e(k-3) = xi(k) + 0.5xi(k-1) + 0.2xi(k-2)


    展开全文
  • MATLAB编程实现产生白噪声与有色噪声序列,代码可读性高,注释到位
  • 资源部含有matlab代码,调用randn函数生成高斯白噪声样本,然后进行统计分析,并且与理论的分布进行了对比,有问题可以私信,有问必答,欢迎你的购买
  • 直接扩频序列调制是用速率很高的伪噪声序列与信息码序列模二相加(波形相乘)后得到复合码序列,用复合码序列去控制载波相位,从而获得直接扩频序列信号的。直接扩频通信具有低截获概率、抗干扰能力强以及易于实现...

    前言

    直接扩频序列调制是用速率很高的伪噪声码序列与信息码序列模二相加(波形相乘)后得到复合码序列,用复合码序列去控制载波相位,从而获得直接扩频序列信号的。直接扩频通信具有低截获概率、抗干扰能力强以及易于实现码分多址等优点,在抗干扰通信及民用移动通信中都得到了广泛的应用。

    文章迭代更新

    文章进行了大改
    发现了为什么振幅是0.2和0.4时,系统无法正常工作的原因。因为bitMultiple函数把运算结果格式转换成了int8,导致精度大量失真!当时转换格式是为内存空间与运行速度做打算,结果今天发现做了负优化!修改后,发现载波振幅对于系统的误码率曲线几乎无影响。
    另外,解调函数也做了修改,会根据输入自动计算判决阈值,解决了之前人为设定阈值
    的局限,阈值采用正态分布解算法,详见demodulate函数内注释。
    新版修改了以下函数:
    1. bitMultiple
    2. demodulate
    3. 新增arrayGroupSum函数

    仿真流程图

    用户1
    扩频
    加扰
    调制
    高斯信道
    walsh码
    扩频
    M序列
    加扰
    载波
    调制
    高斯白噪声
    高斯信道
    用户2
    相加
    解调
    去扰
    用户1解扩
    用户2解扩
    还原信号1
    还原信号2

    关键技术细节

    1. 仿真中,用户码元使用双极性(1和-1)码。
    2. 代码的最后通过还原信号与码元做对比,计算误码率,来评价通信质量的好坏。
    3. 为了研究信噪比对误码率的影响,代码中用了多线程,针对不同的信噪比进行解调去扰解扩,计算不同信噪比下的误码率。
    4. 为了研究载波振幅对信噪比——误码率曲线的影响,增加了振幅的迭代。

    扩频解扰解扩部分

    1. 采用64阶的walsh码矩阵,由walsh函数产生。
    2. 两个用户使用同一扩频矩阵的不同相位(实际上是不同一行)进行扩频。
    3. 扩频时,首先对用户码元按照扩频增益,进行周期延拓,然后点乘扩频码,实现扩频。
    4. 解扩时,首先对输入的去扰后的码按扩频时一样的相位对应点乘扩频码,然后进行判决。

    加扰去扰部分

    1. 加扰采用的是5阶的M序列,反馈系数为67(八进制)。
    2. 加扰同样采用点乘,输入码元1:1地点乘加扰码。
    3. 因为加扰码具有很高的自相关性,所以去扰时只需执行加扰时一样的操作即可。

    加扰原理图

    加扰原理图

    调制解调部分

    1. 调制采用BPSK(二进制移相键控)调制。
    2. 调制所使用的载波为正弦波,通过对其等间隔采样形成离散样值。
    3. 调制时,一个加扰后的码元对应相乘一串正弦波的周期采样值,把相乘结果汇总后就是调制后的结果
    4. 解调时,把一个码元所对应的一串正弦波的周期采样值与载波采样值进行点乘,根据点乘结果中,正数、负数还是零多,判决这是什么码元。

    调制原理图

    调制原理图

    解调原理图

    解调原理图

    高斯信道部分

    1. 使用matlab自带的函数awgn(input,snr(dB),inputPower)
    2. 只要传入需要加高斯噪声的信号,信噪比还有输入信号的功率,即可返回添加了高斯白噪声的信号。

    实验结果

    扩频增益为10时,walsh矩阵为64阶时,不同振幅下信噪比误码率曲线

    扩频增益为10时,walsh矩阵为64阶时,不同振幅下信噪比误码率曲线

    扩频增益为20时,walsh矩阵为64阶时,不同振幅下信噪比误码率曲线

    扩频增益为20时,walsh矩阵为64阶时,不同振幅下信噪比误码率曲线
    以上两次实验结果显示,载波振幅对于系统性能无明显影响,这是文章改版前的一个很大的错误!在此笔者向各位读者道歉。

    误码率随信噪比的变化曲线

    误码率随信噪比的变化曲线
    观察以上曲线可以发现,随着信噪比的上升,误码率以近似于抛物线的形式在不断下降,并最终等于0,这足以证明本次直接序列扩频系统仿真符合预期结果。

    源代码

    可以修改的参数

    1. 轮回次数。通过修改轮回次数,获得更好的曲线,或者更快的代码运行速度,代码运行参考速度:i7-4790:实测速率9.8秒/轮回 @1280*2&1280码元20&10扩频增益。
    2. 用户1的码元数量与扩频增益、用户2的码元数量与扩频增益。修改时注意两用户的码元必须是walsh码阶数的增数倍,同时要保证:用户1的码元数量 * 用户1的扩频增益=用户2的码元数量 * 用户2的扩频增益。
    3. walsh码的阶数。但walsh码阶数必须大于扩频增益。
    4. 用户1扩频相位、用户2扩频相位。注意相位必须1到64之间取值,两个用户的取值不能一样。
    5. M序列阶数与返回系数。反馈系数取值有相关要求,详见百度,作为参数传入时要以八进制表示。
    6. 半个载波周期的采样点数。
    7. 最小最大信噪比,以及其尝试步进。
    8. 最小最大载波振幅,以及其尝试步进。

    !下载链接!

    链接:https://pan.baidu.com/s/14KNoWlD7UyC9WAEza63W5A
    提取码:qg1y

    核心代码

    主函数

    main.m

    %{
        本函数是整个仿真的主函数,用于研究在移动通信时,信噪比对误码率的影响
    首先生成随机双极性码,然后经过扩频,加扰,BPSK调制,加高斯白噪声,混合
    然后模拟接收端的解调,去扰,解扩,判决。
        通过比较接收端判决输出与原来的码元,计算出误码率,首先通信系统的仿真
    以及误码率-信噪比的变化曲线的绘制.
    2019/11/22(以上)
        在当前版本中,添加了对信噪比&误码率曲线随BPSK调制载波振幅变化而变化
    的相关研究功能。因为实际发现,之前的代码中曲线会在-20dB到0dB处出现平缓现象,
    实际排查发现是载波的振幅导致的。振幅研究结果发现:振幅小于0.4时,系统无法正常
    工作!其误码率曲线在后期随着信噪比的增加呈现反常膨胀!
    ^-^想要得到抛物线,只有在载波振幅为1的时候。^-^,不然随着振幅的增加,曲线曲线趋于平缓的信噪比范围会越大;
    同时发现大于0.6的曲线会近似相交于某一点,交点前同等信噪比下振幅高的误码率低,但在交点后
    ***同等的信噪比下振幅小的误码率反而低!!!***
        在本版本中,walsh矩阵的所有码元都被用于扩频,没有使用前一版的矩阵截取方法,
    规避了鬼魅版的非严格正交问题,但前面的截取版本恰恰说明非严格正交对本系统影响曲线
    的影响并不是很大,但却可以使得系统能在更低的误码率下工作。
        至此,代码已经经历10次版本迭代。
    2019/12/1(以上)
        第二次答辩,仍出现部分平滑问题,并且老师说系统性能太好了,与实际系统不符,并且
    曲线不应该受振幅的影响,基于此再对代码做修改,计算了输入信号的平均功率,作为第3
    个参数传给awgn,否则awgn会把输入信号的功率视为0dBW!现在曲线基本重合。但0.20.4
    的问题仍未解决!
        至此,程序第11次版本迭代
    2019/12/2(以上)
        第十二次迭代。
        发现了为什么振幅是0.20.4时,系统无法正常工作的原因,是因为
    bitMultiple函数把运算结果格式转换成了int8,导致精度大量失真!当时转换格式是为
    内存空间与运行速度做打算,结果今天发现做了负优化!修改后,发现载波振幅对于系统
    的误码率曲线几乎无影响。
        另外,解调函数也做了修改,会根据输入自动计算判决阈值,解决了之前人为设定阈值
    的局限,阈值采用正态分布解算法,详见demodulate函数内注释。
    2020/7/1(以上)
        i7-4790:实测速率9.8/轮回 @1280*2&1280码元20&10扩频增益
    %}
    clear variable;
    close all;
    
    mulTimes = 5;  %轮回次数
    
    walshOrder = 64;    %walsh码的阶数,必须大于扩频增益
    
    N1 = 1280*2;   %用户1码元数量
    N2 = 1280;   %用户2码元数量
    
    user1SPgain = 10;  %用户1的扩频增益
    user2SPgain = 20;  %用户2的扩频增益
    
    %记录两个用户的walsh相位,注意相位必须在164之间取值
    %两个用户的取值不能一样
    user1Phase = 2; %用户1的相位
    user2Phase = 16; %用户2的相位
    
    %加扰使用的m序列的参数
    mOrder = 5; %级数5级
    feedBack = 67;%反馈系数67
    
    %调整时,半个周期的采样点数
    samplePiont = 4;
    
    %生成需要使用的walsh码
    walshCode = walsh(walshOrder);
    
    %针对低性能电脑做优化,采取以时间换取空间的思路
    maxSnr = 20;    %尝试的最大信噪比
    minSnr = -20;   %尝试的最小信噪比
    div = 1;      %信噪比的尝试步进
    maxTime = (maxSnr-minSnr)/div;  %尝试次数
    timesUser1Acc = zeros(mulTimes,maxTime);
    timesUser2Acc = zeros(mulTimes,maxTime);
    
    %生成双极性码片
    user1 = genBipolar(N1);
    user2 = genBipolar(N2);   
    
    %扩频
    spread1 = spreadSpectrum(user1,walshCode,user1SPgain,user1Phase);
    spread2 = spreadSpectrum(user2,walshCode,user2SPgain,user2Phase);
    
    %加扰
    Mseq1 = MseqGen(mOrder,feedBack); %用户1加扰用的m序列
    Mseq2 = MseqGen(mOrder,feedBack); %用户2加扰用的m序列
    user1scarm = scarmbling(spread1,Mseq1);
    user2scarm = scarmbling(spread2,Mseq2);
    
    maxAmp = 1.2;    %尝试的最大载波振幅
    minAmp = 0.2;   %尝试的最小载波振幅
    divAmp = 0.2;      %载波振幅的尝试步进
    maxTimesAmp = floor((maxAmp-minAmp)/divAmp);  %振幅尝试次数
    
    ampRecords1 = zeros(maxTimesAmp,maxTime);
    ampRecords2 = zeros(maxTimesAmp,maxTime);
    
    for amp = 1:maxTimesAmp %测试不同的载波振幅下的曲线
        
        fprintf('目前正在第%d个振幅\n',amp);
    
        %调制BPSK
        %生成载波,两用户使用同一个载波
        carrier = (minAmp + divAmp*(amp-1))*sin(0:(pi/samplePiont):(2*pi-2*pi/samplePiont));
        user1modulate = modulate(user1scarm,carrier);
        user2modulate = modulate(user2scarm,carrier);
        
        %计算载波功率
        power = powerCnt(carrier);
        
        for times = 1:mulTimes
    
            fprintf('目前正在第%d个轮回\n',times);
    
            user1Acc = zeros(1,maxTime);
            user2Acc = zeros(1,maxTime);
    
            parfor index = 1:maxTime 
                snr = minSnr + (index-1)*div; %加在发送信号上的高斯噪声的信噪比(dBW)
                            
                %通过高斯信道,添加高斯噪声
                user1send = awgn(user1modulate,snr,power);
                user2send = awgn(user2modulate,snr,power);
    
                %接收方
                receive = user1send + user2send; %收到混合起来的信号
    
                %解调
                demodulateRes = demodulate(receive,carrier);
    
                %去扰
                user1Descarm = deScarmbling(demodulateRes,Mseq1);
                user2Descarm = deScarmbling(demodulateRes,Mseq2);
    
                %解扩
                user1deDS = deSpreadSpectrum(user1Descarm,walshCode,user1SPgain,user1Phase);
                user2deDS = deSpreadSpectrum(user2Descarm,walshCode,user2SPgain,user2Phase);
    
                %计算误码率
                [~,user1Accuracy] = compare(user1,user1deDS);
                [~,user2Accuracy] = compare(user2,user2deDS);
                user1Acc(index) = 1-user1Accuracy;
                user2Acc(index) = 1-user2Accuracy;
            end
            timesUser1Acc(times,:) = user1Acc;
            timesUser2Acc(times,:) = user2Acc;
        end
        %总结统计多次实验的结果
        for i = 1:maxTime
            user1Acc(i) = mean(timesUser1Acc(:,i));
            user2Acc(i) = mean(timesUser2Acc(:,i));
        end
        ampRecords1(amp,:) = user1Acc;
        ampRecords2(amp,:) = user2Acc;
    end
    
    %误码率随信噪比的变化曲线
    figure(1);
    X1 = (minSnr:div:maxSnr-div);
    semilogy(X1,ampRecords1(5,:),'b');
    xlabel('信噪比(dB)');
    ylabel('误码率');
    title('误码率随信噪比的变化曲线');
    hold on;
    semilogy(X1,ampRecords2(5,:),'g');
    legend('用户1(扩频增益10)','用户2(扩频增益20)');
    
    %用户1振幅不同的误码率随信噪比的变化曲线
    figure(2);
    for i = 1:maxTimesAmp
        semilogy(X1,ampRecords1(i,:));
        hold on;
    end
    xlabel('信噪比(dB)');
    ylabel('误码率');
    title('用户1振幅不同的误码率随信噪比的变化曲线');
    legend('0.2','0.4','0.6','0.8','1.0');
    
    %用户2振幅不同的误码率随信噪比的变化曲线
    figure(3);
    for i = 1:maxTimesAmp
        semilogy(X1,ampRecords2(i,:));
        hold on;
    end
    xlabel('信噪比(dB)');
    ylabel('误码率');
    title('用户2振幅不同的误码率随信噪比的变化曲线');
    legend('0.2','0.4','0.6','0.8','1.0');
    

    walsh码生成

    walsh.m

    % 产生 Walsh函数通用函数
    % 参数N表示Walsh函数阶数,N不是2的幂时,通过向无穷大取整使得所得Walsh阶数为2的幂
    function [walsh]=walsh(N)
        M = ceil(log2(N));
        wn = 0;
        for i = 1:M
            w2n = [wn,wn;wn,~wn];
            wn = w2n;
        end
        wn(wn == 0) = -1;
        walsh = int8(wn);
    end
    

    M序列生成

    MseqGen.m

    %M序列发生器
    %order:阶数
    %setBakc:反馈系数
    function res = MseqGen(order,feedBack)
        res = zeros(2^order-1,order);
        feedBack = Oct2Bin(feedBack);
        feedBack = feedBack(2:6);
        temp = rand(1,order);
        temp(temp < 0.5) = 0;
        temp(temp >= 0.5) = 1;
        res(1,:) = temp;
        for i = 2:2^order-1
            newBit = sum(res(i-1,feedBack == 1));
            newBit = mod(newBit,2);
            for j = 1:order-1
                res(i,j+1) = res(i-1,j);
            end
            res(i,1) = newBit;
        end
        res(res == 0) = -1;
    end
    

    扩频模块

    spreadSpectrum.m

    %扩频函数
    %userCode:需要扩频的用户码元
    %PNseq:用于扩频的随机码
    %gain:扩频增益
    %phase:用户扩频码相位
    function res = spreadSpectrum(userCode,PNseq,gain,phase)
        %首先对码元进行周期延拓
        [~,userCode2] = selfCopy(userCode,gain);
        %计算扩频码的行数
        [lineSize,~] = size(PNseq);
        %对扩频码进行重排行序,使初相位位于第一行
        PN = PNseq(phase:lineSize,:);
        PN = [PN;PNseq(1:phase-1,:)];
        res = bitMultiple(userCode2,PN(:)');
    end
    

    解扩模块

    deSpreadSpectrum.m

    %本函数用于解扩
    %userCode:需要解扩的用户码元
    %PNseq:用于扩频的随机码
    %gain:扩频增益
    %phase:用户扩频码相位
    function res = deSpreadSpectrum(userCode,PNseq,gain,phase)
        [lineSize,~] = size(PNseq);
        PN = PNseq(phase:lineSize,:);
        PN = [PN;PNseq(1:phase-1,:)];
        temps = bitMultiple(userCode,PN(:)');
        %这里是matlab的一个小bug,重排序以列作为索引,我的要求是以行
        %作为索引,所以要行列反写再取转置矩阵
        temps = reshape(temps,gain,length(temps)/gain);
        temps = temps';
        %解扩第二步,码元判决
        res = ones(1,length(temps(:))/gain);
        for i = 1:length(temps(:))/gain
            if sum(temps(i,:)) < 0
                res(i) = -1;
            end
        end
    end
    

    加扰模块

    scarmbling.m

    %加扰函数
    %source:被加扰的信号
    %PNCode:用于加扰的扰码
    %1:1的数据加扰
    function res = scarmbling(source,PNCode)
        res = bitMultiple(source,PNCode);
    end
    

    去扰模块

    deScarmbling.m

    %本函数用于对信号做去扰处理,实际操作与加扰完全一致
    %input:需要去扰的信号
    %PNseq:用于去扰的随机序列
    function res = deScarmbling(input,PNseq)
        res = bitMultiple(input,PNseq);
    end
    

    调制模块

    modulate.m

    function res = modulate(source,carrier)
    %调制函数
    %source:被调制信号
    %carrier:载波信号
        sizeSource = length(source);
        sizeCarrier = length(carrier);
        res = zeros(sizeSource,sizeCarrier);
        for i = 1:sizeSource
            res(i,:) = double(source(i))*carrier;
        end
        %把res按行拼接成一行
        res = res';
        res = res(:);
    end
    

    解调模块

    demodulate.m

    function res = demodulate(input,carrier)
    %本函数实现对接收信号的解调
    %原理:极性比较
    %input:被解调的信号
    %carrier:载波信号
    
    %按位循环相乘
    res_bitMultiple = bitMultiple(input,carrier);
    %分组相加
    res_arrayGroupSum = arrayGroupSum(res_bitMultiple,length(carrier));
    %{
        以上处理后,结果正负侧近似于泊松分布,为了处理方便,根据大数定理,
        近似地视为正态分布,从而计算阈值。按照正态分布,-1,1,0平分概率,每一个分得1/3的概率。
    %}
    %计算标准差
    res_std = std(res_arrayGroupSum);
    %计算平均数
    res_mean = mean(res_arrayGroupSum);
    %解算阈值
    syms temp;
    %以下计算时会有一个"Unable to solve symbolically. Returning a numeric solution using
    %vpasolve"的警告,让系统不显示
    warning off;
    threshhold_negative = double(solve(normcdf(temp,res_mean,res_std)==1/3));
    %计算概率为2/3时的阈值,也就是判决为1的阈值
    threshhold_postive = double(solve(normcdf(temp,res_mean,res_std)==2/3));
    res = zeros(1,length(res_arrayGroupSum));
    for i = 1:length(res_arrayGroupSum)
        if res_arrayGroupSum(i) > threshhold_postive
            res(i) = 1;
        elseif res_arrayGroupSum(i) < threshhold_negative
            res(i) = -1;
        end
    end
    end
    

    模块测试代码

    本部分的代码用于单独测试三大模块是否正常工作

    扩频解扩模块测试

    testSpreadSpectrum.m

    %本代码用于测试扩频,解扩模块是否正常工作
    %注意扩频与解扩针对的是双极性码
    function testSpreadSpectrum()
        raw = genBipolar(64);
        walshCode = walsh(64);
        afterDSSS = spreadSpectrum(raw,walshCode,4,2);
        res= deSpreadSpectrum(afterDSSS,walshCode,4,2);
        [trueNum,accuracy] = compare(raw,res);
        fprintf("正确码元数量:%d\n正确率:%f\n",trueNum,accuracy); 
    end
    

    加扰去扰模块测试

    testScarmbling.m

    %本代码用于测试加扰,去扰模块能否正常工作
    function testScarmbling()
        %M序列发生器可以自定义阶数和反馈系数
        Mseq = MseqGen(5,67); %产生加扰用的m序列
        raw = tripleGen(1e6);
        afterScarmb = scarmbling(raw,Mseq);
        res = deScarmbling(afterScarmb,Mseq);
        [trueNum,accuracy] = compare(raw,res);
        fprintf("正确码元数量:%d\n正确率:%f\n",trueNum,accuracy);    
    end
    

    调制解调模块测试

    testModulate.m

    %本代码用于测试BPSK调制与解调模块能否正常工作
    function testModulate()
        carrier = 10*sin(0:pi/32:2*pi-pi/32);
        raw = tripleGen(1e6);
        afterModu = modulate(raw,carrier);
        res = demodulate(afterModu,carrier);
        [trueNum,accuracy] = compare(raw,res);
        fprintf("正确码元数量:%d\n正确率:%f\n",trueNum,accuracy);   
    end
    

    非核心代码

    八进制转二进制模块

    Oct2Bin.m

    %八进制转二进制,返回一个二进制数组
    function res = Oct2Bin(value)
        [bitNum,bit] = getPalces(value);
        res = int8(zeros(1,bitNum*3));
        for i = 1:bitNum*3
            bitMain = bit(floor((i-1)/3)+1);
            switch mod(i-1,3)
                case 0
                    if bitMain >= 4
                        res(i) = 1;
                    end
                case 1
                    if bitMain == 2 || bitMain == 3 || bitMain == 6 || bitMain == 7
                       res(i) = 1;
                    end
                case 2
                    if mod(bitMain,2) ~= 0
                        res(i) = 1;
                    end
            end
        end
    end
    

    点乘模块

    bitMultiple.m

    %按bit循环相乘函数,如果信号长度不一致,则较短的信号被循环使用
    %signal1:相乘信号1
    %signal2:相乘信号2
    %res:相乘结果
    function res = bitMultiple(signal_1,signal_2)
        sizeSource = length(signal_1);
        sizePNCode = length(signal_2);
        %如果长度不满足条件时,调换顺序递归调用
        if sizeSource < sizePNCode
            res = bitMultiple(signal_2,signal_1);
            return;
        end
        res = zeros(1,sizeSource);
        for i = 1:sizeSource
            res(i) = signal_1(i) * signal_2(mod(i-1,sizePNCode)+1);
        end
        res = int8(res);
    end
    

    码元比较模块

    compare.m

    %比较两个数组
    %input1:数组1
    %input2:数组2
    %res:正确码元数量
    %accuracy:正确率
    function [res,accuracy] = compare(input1,input2)
        temp = (input1 == input2);
        res = length(find(temp == 1));
        sizeInput = max(length(input1),length(input2));
        accuracy = res/double(sizeInput);
    end
    

    双极性码生成模块

    genBipolar.m

    %产生双极性码
    %num:双极性码的规模
    %res:满载双极性码的数组
    function res = genBipolar(num)
        res = rand(1,num);
        res = value2Bipolar(res);
    end
    

    数字位数获取模块

    getPlaces.m

    %统计一个数字的位数,并返回每位的数字
    function [Places,item] = getPalces(value)
        Places = 0;
        temp = abs(value);
        while temp ~= 0
            temp = floor(temp/10);
            Places = Places + 1;
        end
        item = uint8(zeros(1,Places));
        temp = value;
        for i = Places:-1:1
            item(i) = mod(temp,10);
            temp = floor(temp/10);
        end
    end
    

    功率计算模块

    powerCnt.m

    %计算输入信号的平均功率(dBW)
    function res = powerCnt(input)
        res = 10*log10(sum(input.*input));
    end
    

    周期延拓模块

    selfCopy.m

    %本函数实现数组的周期延拓
    %input:需要延拓的数组
    %times-1:延拓的次数
    %比如selfCopy([1,2,3],2) = [1,2,3,1,2,3,1,2,3]
    function [res,res2] = selfCopy(input,times)
        if times <= 1
            res = input;
        else
            res = zeros(length(input),times);
            for i = 1:times
                res(:,i) = input;
            end
            res2 = res';
            res2 = res2(:)';
            res = res(:)';
        end
    end
    

    三极码生成模块

    tripleGen.m

    %本函数用于生成指定数量的随机三极性码,主要用于做模块测试
    %num:数量
    %res:满载结果的数组
    function res = tripleGen(num)
        raw = rand(1,num);
        for i = 1:num
            if raw(i) >= 0.75
                raw(i) = 1;
            elseif raw(i) <= 0.25
                raw(i) = -1;
            else
                raw(i) = 0;
            end
        end
        res = int8(raw);
    end
    

    0~1转双极性模块

    value2Bipolar.m

    % 本函数用户把随机生成的0~1之间的double数转换成双极性码
    function res = value2Bipolar(input)
        res = zeros(1,length(input));
        for index = 1:length(input)
            if(input(index) < 0.5)
                res(index) = -1;
            else
                res(index) = 1;
            end
        end
        res = int8(res);
    end
    

    分组相加模块

    arrayGroupSum.m

    function res = arrayGroupSum(input,group_num)
    %arrayGroupSum 数组分组求和,如果数组长度除去group_num无法整除,则截短input数组
    %input 目标数组
    %group_num 分组数
    len = floor(length(input)/group_num);
    res = zeros(1,len);
    for i = 1:len
        res(i) = sum(input((i-1)*group_num+1:i*group_num));
    end
    end
    
    展开全文
  • matlab产生高斯噪声

    万次阅读 2018-03-09 19:52:59
    randn randn(random normal distribution)是一种产生标准正态分布的随机数或矩阵的函数,属于MATLAB函数。返回一个n*n的随机项的矩阵。如果n不是个数量,将返回错误信息。MATLAB函数randn简介编辑用法:Y = ...

    randn

     
    randn(random normal distribution)是一种产生标准 正态分布随机数矩阵的函数,属于MATLAB函数。返回一个n*n的随机项的矩阵。如果n不是个数量,将返回错误信息。

    MATLAB函数randn简介

    编辑
    用法:
    Y =  randn(n)
    返回一个n*n的随机项的矩阵。如果n不是个数量,将返回错误信息。
    Y = randn(m,n) 或 Y = randn([m n])
    返回一个m*n的随机项矩阵。
    Y = randn(m,n,p,...) 或 Y = randn([m n p...])
    产生随机 数组
    Y = randn(size(A))
    返回一个和A有同样 维数大小的 随机数组。
    randn
    返回一个每次都变化的数量。


    应用举例


    Example 1. R =  randn(3,4) 将生成 矩阵
    R =
    1.1650 0.3516 0.0591 0.8717
    0.6268 -0.6965 1.7971 -1.4462
    0.0751 1.6961 0.2641 -0.7012
    For a histogram of the randn distribution, see hist.
    Example 2. 产生一个随机分布的指定 均值方差的矩阵:将randn产生的结果乘以 标准差,然后加上期望均值即可。例如,产生均值为0.6,方差为0.1的一个5*5的 随机数方式如下:
    x = .6 + sqrt(0.1) *  randn(5)
    x =
    0.8713 0.4735 0.8114 0.0927 0.7672
    0.9966 0.8182 0.9766 0.6814 0.6694
    0.0960 0.8579 0.2197 0.2659 0.3085
    0.1443 0.8251 0.5937 1.0475 -0.0864
    0.7806 1.0080 0.5504 0.3454 0.5813
    Example 3.设置randn到其默认的初始状态格式:randn('state', 0);每次初始化randn到不同的状态的格式randn('state', sum(100*clock))。
    如何用randn产生两个相同的矩阵:
    s = randn('state');%保存当前状态,
    u1 = randn(100);%产生100个值,
    randn('state',s);%复位状态,
    u2 = randn(100); %并重复序列。格式调用如下:
    s = randn('state');
    u1 = randn(3)
    randn('state',s);
    u2 = randn(3)
    运行结果如下:
    u1 =1.6039 0.8957 -0.1072
    -0.3728 -1.3913 -0.8276
    0.7578 0.3116 0.0592
    u2 =1.6039 0.8957 -0.1072
    -0.3728 -1.3913 -0.8276
    0.7578 0.3116 0.0592
    从结果可以看出,这两个矩阵的各个元素均相同,所以这两个就是等价矩阵。
    其他类似函数: rand, randperm, sprand, sprandn
    Y = randn(n)
    Y = randn(m,n)
    Y = randn([m n])
    Y = randn(m,n,p,...)
    Y = randn([m n p...])
    Y = randn(size(A))
    randn
    s = randn('state')
    展开全文
  • 实验一 基于 MATLAB 的白噪声信号 u(n)sinc 函数chirp 信号产生实验 一 实验目的 1. 学会使用 MATLAB 2. 通过实验了解 MATLAB 如何产生各种常用信号 3. 掌握 MATLAB 的编程方法 二 实验内容 1. 用 MATLAB 编程产生一...
  • matlab产生高斯白噪声

    千次阅读 多人点赞 2021-02-06 17:57:32
    (2) randn:产生均值为0、方差为1的高斯白噪声。 (3) randperm(n):产生1到n的均匀分布随机序列。 (4) normrnd(a,b,c,d):产生均值为a、方差为b大小为cXd的 随机矩阵。 rand:返回一个在区间 (0,1) 内均匀...

    函数介绍

    matlab里和随机数有关的函数:
    (1) rand:产生均值为0.5、幅度在0~1之间的伪随机数。
    (2) randn:产生均值为0、方差为1的高斯白噪声。
    (3) randperm(n):产生1到n的均匀分布随机序列。
    (4) normrnd(a,b,c,d):产生均值为a、方差为b大小为cXd的 随机矩阵。

    • rand:返回一个在区间 (0,1) 内均匀分布的随机数。

      • rand(n):生成0到1之间的n阶( n×n )随机数方阵。
      • rand(m,n):生成0到1之间的m×n的随机数矩阵。
    • randn:返回一个从标准正态分布中得到的随机标量。

      • randn()命令是产生白噪声的,白噪声应该是0均值,方差为1的一组数。
      • 同rand函数一样,randn(n),randn(m,n)含义与上述一致。
      • randn(size(A)),返回一个和A有同样维数大小的随机数组。
    • randperm:整数的随机排列。

      • p = randperm(n) 返回行向量,其中包含从 1 到 n 没有重复元素的整数随机排列。
      • p = randperm(n,k) 返回行向量,其中包含在 1 到 n 之间随机选择的 k 个唯一整数
    • normrnd:生成服从正态分布的随机数

      • r = normrnd(mu,sigma)均值参数为 mu标准差参数为 sigma 的正态分布中生成随机数。
      • R=norrmrnd(MU,SIGMA,m):从均值参数为 mu标准差参数为 sigma 的正态分布中生成随机数,矩阵的形式由m定义。m是一个1×2向量,其中的两个元素分别代表返回值R 中行与列的维数。
      • R=normrnd(MU,SIGMA,m,n): 生成m×n形式的正态分布的随机数矩阵。

    注:

    • 一般来说,可以使用公式r = a + (b-a).*rand(N,1)生成区间 (a,b) 内的 N 个随机数。
    • rand是0-1的均匀分布,randn是均值为0方差为1的正态分布。
    • 理论上randn()生成的随机数分布范围为(-∞,+∞),即无穷大。Matlab中randn()是产生正态分布的随机数或矩阵的函数,它产生均值为0,方差为1,标准差为1的正态分布的随机数或矩阵的函数。

    高斯白噪声函数

    • 高斯白噪声概念解释
      • 高斯白噪声(white Gaussian noise; WGN):均匀分布于给定频带上的高斯噪声
        • 如果一个噪声,它的幅度服从高斯分布,而它的功率谱密度又是均匀分布的,则称它为高斯白噪声。
        • 高斯白噪声中的高斯是指概率分布是正态函数,而白噪声是指:它的二阶矩不相关,一阶矩为常数,是指先后信号在时间上的相关性。这是考察一个信号的两个不同方面的问题。
        • 热噪声散粒噪声高斯白噪声
    • matlab高斯白噪声函数介绍:——wgn( )、awgn( )
      • WGN:产生高斯白噪声
        • y = wgn(m,n,p) 产生一个m行n列的高斯白噪声的矩阵,p以dBW为单位指定输出噪声的强度
        • y = wgn(m,n,p,imp)欧姆(Ohm)为单位指定负载阻抗
        • y= **wgn(…,outputtype)**指定输出类型,OUTPUTTYPE可以是’real’或’complex’来获得复噪声信号。
        • y=**wgn(…,powertype)**指定p的单位。POWERTYPE可以是’dBW’, ‘dBm’或’linear’。线性强度(linear power)以瓦特(Watt)为单位。
        • linear表示线性强度(linear power),单位为Watt。如果输入其他:‘dBw’或缺省则表示用dBw作为功率单位
      • AWGN:在某一信号中加入高斯白噪声
        • y = awgn(x,SNR) 在信号x中加入高斯白噪声。信噪比SNR以dB为单位。x的强度假定为0dBW。如果x是 复数,就加入复噪声。
        • y = awgn(x,SNR,SIGPOWER) 如果SIGPOWER是数值,则其代表以dBW为单位的信号强度;如果SIGPOWER为’measured’,则函数将在加入噪声之前测定信号强度。
        • y = awgn(x,SNR,SIGPOWER,STATE) 重置RANDN的状态。
        • y = awgn(…,POWERTYPE) 指定SNR和SIGPOWER的单位。POWERTYPE可以是’dB’或’linear’。如果POWERTYPE是’dB’,那么SNR以dB为单位,而SIGPOWER以dBW为单位。如果POWERTYPE是’linear’,那么SNR作为比值来度量,而SIGPOWER以瓦特为单位。

    注:

    • 一阶矩就是随机变量的期望二阶矩就是随机变量平方的期望,以此可以类推高阶的矩。
    • 分贝(decibel, dB):分贝(dB)是表示相对功率或幅度电平的标准单位,换句话说,就是我们用来表示两个能量之间的差别的一种表示单位,它不是一个绝对单位。
      • 例如,电子系统中将电压、电流、功率等物理量的强弱通称为电平,电平的单位通常就以分贝表示,即事先取一个电压或电流作为参考值(0dB),用待表示的量与参考值之比取对数(以10为底的对数),再乘以20作为电平的分贝数(功率的电平值改乘10)。
      • 分贝瓦(dBW, dB Watt):指以1W的输出功率为基准时,用分贝来测量的功率放大器的功率值
        • dBm (dB-milliWatt):即与1milliWatt(毫瓦)作比较得出的数字。
        • 0 dBm = 1 mW 10 dBm = 10mW 20 dBm = 100 mW

    总结

    • 在matlab中无论是wgn还是awgn函数实质都是由randn函数产生的噪声
    • 即:wgn函数中调用了randn函数,而awgn函数中调用了wgn函数。
    • awgn(x,snr,’measured’,'linear’)表示向已知信号添加某个信噪比(SNR)的高斯白噪声“,命令的作用是对原信号x添加信噪比(比值)为SNR的噪声,在添加之前先估计信号x的强度
      • 定义解释:
        • SNR就是信号的强度除以噪声的强度或者信号功率与噪声功率之比((注:由于采用的是比值而非db,所以与下面“计算信噪比”所使用的方式不同,即没有求对数步骤))
        • 信号的强度指的就是信号的能量,在连续的情形就是对x平方后求积分,而在离散的情形是求和代替积分(在matlab中sigPower= sum(abs(sig(: )).^2)/length(sig(: )),这就是信号的强度,这里sig(: )为信号。)
      • 在求出x的强度后,结合指定的信噪比,就可以求出需要添加的噪声的强度noisePower = sigPower/ SNR。由于使用的是高斯白噪声即randn函数,而randn的结果是一个强度为1的随机序列(自己试试sum(randn(1000,1).^2)/1000就知道了,注意信号的长度不能太小)。于是,所要添加的噪声信号显然就是: sqrt(noisePower)*randn(n,1)其中n为信号长度。
      • matlab 例子:
    X = sqrt(2)*sin(0:pi/1000000:6*pi);               % 产生正弦信号
    Y = awgn(X,10,'measured');                        % 加入信噪比为10db的噪声,加入前预估信号的功率(强度)
    sigPower = sum(abs(X).^2)/length(X)           	  % 求出信号功率
    noisePower=sum(abs(Y-X).^2)/length(Y-X)  		  % 求出噪声功率
    SNR=10*log10(sigPower/noisePower)        		  % 由信噪比定义求出信噪比,单位为db
    
    • 运算结果如下所示:
      在这里插入图片描述
      • WGN(m,n,p)产生功率为p dBW的m*n的高斯白噪声矩阵,其中p是以dbW为单位的输出强度。若要产生一个均值0,方差为0.0965 的高斯白噪声,不可直接用WGN(N,1,0.0965)产生(单位不对应)
      • 对高斯白噪声,其方差和功率(单位为W)是一样的。因此,对方差,要做的只是将w变换成dbw,即dbw=10log(w)
      • 做法如下有两种:
    %% 方法一:
    N=1000;
    x=sqrt(0.0965)*randn(N,1);
    Px=(x.'*x)/N   % 验证,这里Px的求法与上面noisePower=sum(abs(Y-X).^2)/length(Y-X)的求法是一致的
    %% 方法二:
    N=1000;
    y=wgn(N,1,10*log10(0.0965));
    Py=(y.'*y)/N   % 验证
    

    信噪比,英文名称叫做SNR或S/N(Signal Noise Ratio),是指系统中信号与噪声的比例。信号指的是来自设备外部需要通过这台设备进行处理的电子信号,噪声是指经过该设备后产生的原信号中并不存在的无规则的额外信号(或信息),并且该种信号并不随原信号的变化而变化。
    信噪比的计量单位是dB,其计算方法是10LOG(Ps/Pn),其中Ps和Pn分别代表信号和噪声的有效功率,也可以换算成电压幅值的比率关系:20LOG(Vs/Vn),Vs和Vn分别代表信号和噪声电压的“有效值”。信噪比应该越高越好。

    参考来源

    matlab 中产生高斯白噪声
    高斯白噪声及Matlab常用实现方法
    关于dB 分贝
    Matlab产生高斯白噪声
    MATLAB产生特定功率谱密度的高斯白噪声的两种方法

    展开全文
  • 这里有三段程序,分别是产生高斯白噪声的程序,信号加载高斯白噪声的程序,产生有色噪声的程序。是本人搜集的,特此分享。
  • 本文主要使用MATLAB产生一均匀分布的白噪声信号u(n),画出其波形,并画出其直方图,检验其分布情况。 rand函数默认产生的是均值为0.5,幅度在0~1之间均有分布的伪随机数。 代码如下: %产生一均匀分布的白噪声...
  • MATLAB产生高斯白噪声

    热门讨论 2009-12-06 20:53:01
    wgn 高斯白噪声 MATLAB 这个经过修改的函数,
  • % p 如果比 0.05 小就不是白噪声序列,可以使用时间序列 % 某种现象的指标数值按照时间顺序排列而成的数值序列 % 因为时间序列是某个指标数值长期变化的数值表现 % 所以时间序列数值变化背后必然蕴含着数值变换的规
  • Matlab时间序列分析

    万次阅读 多人点赞 2018-11-13 18:53:46
    在引入时间序列前,先介绍几个matlab函数 matlab中的gallery函数简析 Matlab 中的 gallery 函数是一个测试矩阵生成函数。当我们需要对某些算法进行测试的时候,可以利用gallery函数来生成各种性质的测试矩阵。其用法...
  • Matlab产生高斯白噪声

    千次阅读 2017-11-03 22:55:42
     在matlab中无论是wgn还是awgn函数,实质都是由randn函数产生噪声。即:wgn函数中调用了randn函数,而awgn函数中调用了wgn函数。  根据awgn的实现代码可以知道”向已知信号添加某个信噪比(SNR)的高斯白噪声...
  • matlab生成m序列的方法

    千次阅读 2020-05-31 10:14:01
    文章目录matlab生成m序列的方法1.m序列基本知识点2.matlab产生m序列2.1根据产生原理编写生成函数2.1.1生成m序列的函数:2.1.2调用已编写函数生成m序列2.2利用idinputidinputidinput函数 引言 m序列属于伪随机序列的...
  • MATLAB中白噪声产生

    千次阅读 2015-09-23 10:02:00
    rand产生的是[0,1]上的均匀分布的随机序列randn产生均值为0,方差为1的高斯随机序列,也就是白噪声序列 rand产生的是均匀分布白噪声序列randn产生的是正态分布的白噪声序列 MATLAB还提供了两个产生高斯白噪声的...
  • 它的灵感来自于 matlab 函数 fitdist。 ----------- ----------- 这是脚本的第一个版本,因此,很快就会有一些变化。 我没有进行任何新的工作。 所有的功劳都归于 [1] 和 [2]。 任何改进脚本的评论或提议都受到热烈...
  • 在C环境下,实现高斯白噪声序列产生,使用瑞利分布
  • matlab产生制定均值和方差的白噪声

    千次阅读 2019-08-24 10:14:19
    w =a+ b.*randn(m,n); 其中:a为均值; b为标准差; m为需要产生几行; n为需要产生几列. 验证方式为:mean(w) var(w) ...
  • MATLAB中高斯白噪声产生

    千次阅读 2020-06-15 10:46:28
    Number = 100000; noise = 1/sqrt(2)*randn(1,...复高斯白噪声产生 randn(1,Number)产生Number点均值为0,方差为1的高斯白噪声 代码中的1/sqrt(2)使得复高斯白噪声的方差(即功率为1)。 数学上看,noise=N(0,.
  • 噪声产生MATLAB实现

    千次阅读 2021-01-03 10:35:01
    噪声
  • 该函数最初是为去除时间序列水速数据中的尖峰噪声而编写的,但可用于通用目的。 基本思想来自 Goring 和 Nikora (2002),它考虑了时间序列信号的一阶和二阶导数。 请参阅参考资料中的详细信息。 参考- Mori, N.、T....
  • 这是我在无线通信学习中做的一个Matlab程序,流程是:先生成一个随机序列作为信号,然后产生了一个长度为15的m序列用作用户1的扩频序列,循环移位作为用户2的扩频序列。经过模拟载波调制后,加入0dB信噪比的高斯白...
  • 噪声和M序列的生成

    2018-04-08 17:02:14
    高斯白噪声和M序列的生成,用于系统辨识的输入信号,MATLAB编程实现
  • matlab产生高斯白噪声

    万次阅读 2017-10-24 10:37:47
    (2) randn:产生均值为0、方差为1的高斯白噪声 (3) randperm(n):产生1到n的均匀分布随机序列 (4) normrnd(a,b,c,d):产生均值为a、方差为b大小为cXd的随机矩阵 rand rand(n):生成0到1之间的n阶随机数方阵...
  • Matlab - 产生高斯噪声

    万次阅读 2015-10-05 14:03:07
    % MATLAB 命令是normrnd。 %1)R=normrnd(MU,SIGMA):生成服从正态分布(MU参数代表均值,DELTA参数代表标准差)的随机数。 % 输入的向量或矩阵MU和SIGMA必须形式相同,输出R也和它们形式相同。
  • MATLAB产生高斯白噪声

    千次阅读 2019-06-13 10:06:21
    2) randn产生均值为0,方差为1的高斯随机序列,也就是白噪声序列; 也就是说,可以直接使用上面两个函数对原始信号添加噪声(例如y=x+rand(length(x),1)或者y=x+randn(length(x),1)) ============================...
  • Matlab噪声高斯噪声

    2020-09-18 15:04:20
    实现书本《随机控制》上关于生成高斯白噪声的方法。 白噪声就是标准均匀分布伪随机数列。 1.标准均匀分布函数,均值1/2,方差1/12; x1=1973; y=zeros(1,500); for i=1:500 x1=mod(91*x1,10^4); y(1,i)=x1/10000...
  • PRBS 基于 3 到 9 位内存延迟并输出 2^3 -1 到 2^9 -1 二进制序列长度。 用途:这些发生器可用于过程识别作为命令信号的噪声源。
  • % 离散序列的功率谱% 具有噪音检验的一维序列x(n)的离散功率谱分析,ol(lw)频率,tl(lw)周期,sl(lw)离散功率谱,% st95(lw)红噪音或白噪音谱的95%置信上限,其中lw=[n/2.]。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,672
精华内容 1,868
关键字:

matlab产生噪声序列

matlab 订阅