精华内容
下载资源
问答
  • 2019-11-10 19:50:37

    空间平滑MUSIC算法(只适用于均匀线性阵列)

    1、由于信号阵列会接收到不同方向上的相干信号,而相干信号会导致信源协方差矩阵的秩亏缺(不是满秩矩阵),从而使得信号特征向量发散到噪声子空间去(就是组成噪声子空间的向量中混入了信号源对应的特征向量)。

    2、对于相干信号波达方向估计:就是想办法把信号协方差矩阵的秩恢复到等于信号源的个数(这样我们得到的信号子空间的秩才是和信号源的个数相等的)。

    3、其中一种思路:在进行空间阵列谱函数估计出波达方向之前,先进行预处理,将协方差的秩恢复到信号源数(即去相关)

    4、去相关中的一种思路:降维处理,比如空间平滑MUSIC算法,牺牲有效阵列孔径来实现信号源的去相干(因为毕竟子阵列的孔径比原阵列的孔径小)

                                            

    clear all;close all;clc
    %阵列的基础信息
    derad=pi/180;%角度转换成弧度
    redeg=180/pi;%弧度转换成角度
    twpi=2*pi;
    Melm=7;%阵元数目
    K=6;%阵元数目(后面会从这两种数目中选择一个)
    dd=0.5;%阵元间距
    d=0:dd:(Melm-1)*dd;%阵元总长度(7-1)*0.5=3,这里共7个元素!
    iwave=3;%信号源数目(有用信号,不包括噪声源)
    theta=[0 30 60];%这三个信号源的入射角度分别是0、30、60
    n=200;%采样点数目(即快拍数)
    A=exp(-j*twpi*d.'*sin(theta*derad));%构造方向矢量(阵列流型),注意d.’
     
    %构造相干信号源(了解相干信号的定义)
    S0=randn(iwave-1,n);%信号源符合高斯分布(正态分布),矩阵维度自行设定,主要考虑后面运算不要起冲突就行
    S=[S0(1,:);S0];%S0是2*n维,S0(1,:)是1*n维,而;是下一行,所以构成了相干信号源S是3*n维
    X0=A*S;%阵列接收信号(与信号源和方向向量有关)
    snr=10;%信噪比,也可以设定为好几个值,snr=[10 20 30],但要保证后面矩阵运算维数不起冲突
    X=awgn(X0,snr,'measured');%往接收信号中添加噪声,模拟真实环境(存在噪声)接收到的信号
    Rxxm=X*X'/n;%协方差矩阵(这只是待定的,最后要用的是Rxx)
    issp=1;%选用的平滑算法模式,我们设计了几种,从中选用一种,比如issp=1是第一种模式,issp=2是第二种
     
    %设计了几种空间平滑算法,以子程序中自定义函数的形式进行调用(ssp和issp)
    if issp == 1%注意,等于是==,而赋值是=
      Rxx = ssp(Rxxm,K);%自定义函数ssp的输出,经过在ssp函数中进行空间平滑算法,输出协方差矩阵
    elseif issp == 2%选择第二种空间平滑算法,注意是在另一个自定义函数issp中运算
      Rxx = issp(Rxxm,K);%将参数Rxxm和kelm代入自定义函数中,最后输出Rxx,比以前那些算法的要好
    else
    Rxx=Rxxm;%没有选用设计的空间平滑算法,所以沿用以前方法的协方差矩阵
    K=Melm;%因为沿用以前方法,所以阵元数目用7
    end
     
    %对协方差矩阵进行特征值分解
    [EV,D]=eig(Rxx);%分解后得到特征值组成的向量形式D和特征值对应特征向量组成的矩阵EV
    DA=diag(D)';%构造对角阵DA,特征值向量D的元素是对角线上的数值,并将矩阵转置(下下句有用)
    [DA,I]=sort(DA);%按照特征值大小对对角阵元素升序排列,原先的顺序为I
    DA=fliplr(DA);%翻转矩阵(中心对称旋转),所以对角阵元素此时是按照降序排列的(结合diag(D)’
    EV=fliplr(EV(:,I));%先使EV中列向量按照I的顺序排列,在进行翻转
     
    %阵列进行扫描
    for jiaodu =1:361%从1循环到361(而且循环数必须为正,>0)
    angle(jiaodu )=(jiaodu -181)/2;%使得阵列扫描角度是从-90到90
    phim=derad*angle(jiaodu );%提前精简,避免下面表达式过于臃肿
    a=exp(-j*twpi*d(1:K)*sin(phim)).';%方向向量,其中取d的前kelm个元素,不会取到Melm个元素。并且要转置.主要是后面矩阵运算,维度有要求
    L=iwave;%信源数目
    En=EV(:,L+1:K);%噪声空间,前L列为信号空间,而后面的则为噪声空间
    P(jiaodu )=(a'*a)/(a'*En*En'*a);%空间阵列谱函数估计
    end
    P=abs(P);%求绝对值或者复数实部与虚部的平方和的算术平方根
    Pmax=max(P);
    P=10*log10(P/Pmax);%归一化,结合信噪比的相关公式
    figure%创建一个用来显示图形输出的一个窗口对象
    h=plot(angle,P);
    set(h,'Linewidth',2)%将所有线宽设置为2
    xlabel('angle(degree)')
    ylabel('magnitude(dB)')
    axis([-90 90 -60 0])%设置坐标轴范围
    set(gca,'XTick',[-90:30:90],'YTick',[-60:10:0])%设置网格的显示格式,gca获取当前figure的句柄,选择x、y轴的要进行标注的刻度
    grid on 
    hold on
    legend('平滑MUSIC')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %子程序1
    %自定义函数ssp
    function Rxx=ssp(Rxxm,K)%自定义函数(参数1待定协方差矩阵,参数2阵元数目)
    [M,MM]=size(Rxxm);%返回Rxxm矩阵的维度大小
    N=M-K+1;%子阵列的数目
    Rxx=zeros(K,K);%初始化矩阵Rxx,类似于初始化为0
    for i=1:N%第1个子阵列到第N个子阵列
    Rxx=Rxx+Rxxm(i:i+K-1,i:i+K-1);%向前滑动,注意每个子阵列是k个阵元,不是k-1个
    end
    Rxx=Rxx/N;%子阵列协方差矩阵可以相加后平均取代原来意义上的协方差矩阵
    %子程序2
    %自定义函数issp
    function Rxx=issp(Rxxm,K)
    [M,MM]=size(Rxxm);
    N=M-K+1;
    J=fliplr(eye(M));%翻转对角阵,中心对称,维数为M*M
    Rxxb=(Rxxm+J*Rxxm.'*J)/2;%前后向平滑协方差矩阵(同时结合了Toeplitz的协方差矩阵估计值);共轭倒序得到的信号子空间与原来的一样
    Rxx=zeros(K,K);
    for i=1:N
    Rxx=Rxx+Rxxb(i:i+K-1,i:i+K-1);
    end
    Rxx=Rxx/N;

    将MUSIC算法的步骤归纳为:

    1、根据n个接收信号矢量,得到协方差矩阵的估计值Rxxm=X*X’/n

    2、若各子阵列的阵列流形相同(等距线阵满足这个假设),则子阵列协方差矩阵可以相加后平均取代原来意义上的协方差矩阵Rxx=Rxx/N

    3、对平均后得到的协方差矩阵进行特征值分解[EV,D]=eig(Rxx)

    4、按照特征值的大小顺序,把与信号个数kelm相等的最大特征值对应的特征向量看作信号子空间,剩下的(M-K)个特征值对应特征向量看作噪声子空间

    5、使θ变化,即扫描角度从-90到90,按照空间阵列谱函数公式来计算,通过寻求峰值来得到波达方向的估计值

     

    更多相关内容
  • 【3】前向平滑技术+MUSIC; 【4】后向平滑技术+MUSIC; 【5】前后双向平滑技术+MUSIC; 【适应对象】: 雷达专业、阵列信号处理专业学生; 【资源特点】: 编程规范,注释明细; 【使用建议】: 此资源为较基础的...
  • 近年来空间平滑技术在减小有限次采样引起的误差、相关信号解等相关方面取得了巨大的进步,此算法与采样协方差矩阵求逆(DMI)法相比,可使信干噪比有一定提高,但改善程度有限。针对空间平滑技术这一特点提出了一种...
  • 改进空间平滑技术的ESPRIT探地雷达时延估计
  • 针对移动通信环境中多径传播的信号到来方向估计问题, 结合空间平滑技术, 提出一种空间平滑的 ESPR IT 算法, 并给出了多径信号数目的推定方法. 仿真结果证明, 与ESPR IT 算法相比, 所提出的算法不仅适用于 ...
  • 系统分析了搜索空间逐步平滑策略的可行性。在离散的情况下,分析了搜索空间逐步平滑法在离散问题中的应用,并与没有经过平滑技术处理过的遗传算法的运行过程进行了比较,得出搜索空间逐步平滑法收敛效率更高。
  • 针对基于阵列协方差矩阵特征分解的子空间类算法存在的问题,提出了一种基于改进空间平滑的新方法。首先介绍了“等效信源”的概念,在此基础上分析了当目标数多于发射阵元数时,一些基于子空间类算法失效的原因;从理论上...
  • 摘 要:讨论了MUSIC空间谱估计算法的改进方法-空间平滑法、超分辨率在米波雷达中的应用,主要是通过建立阵列信号模型将空间平滑法应用于米波雷达实现超角分辨,仿真结果表明空间平滑法可以有效地分辨常规MUSIC法...
  • 为解决在冲击噪声背景和相干信源条件下,高斯白噪声的相干信源DOA(Direction习以rrival)估计算法失效的问题,提出了基于虚拟空间平滑共变系数矩阵(ROC-VSS:RObust Covariation-based virtual spatial Smoothing)的DOA...
  • 首先分析了直接平滑原物理阵列和四阶累积量产生的虚拟阵列不能解相干的原因;然后推导了基于累量域聚焦技术解相干的原理;最后提出了直接累量域聚焦、前后向累量域聚焦和二次前后向累量域聚焦算法,并对其进行了性能...
  • 空间平滑doa估计

    2013-04-05 20:31:54
    空间平滑技术的DOA估计,可以去相关,解决了相干信号问题下DOA 的估计问题
  • 为了能够提升分解矩阵的稀疏表达能力, 提出了一种新的基于平滑l0范数的正交子空间非负矩阵分解方法。通过将分解矩阵的正交性及平滑l0范数约束同时引入矩阵分解的目标函数中一起进行优化, 大大降低了计算复杂度, 并...
  • 研究论文-基于协方差矩阵的空间平滑解相干算法
  • 在充分调研大型激光系统多种光束平滑技术基础上,针对高功率准分子激光系统波长短、宽频带等特点,选用了无阶梯诱导空间非相干(EFISI)技术来实现角多路高功率准分子激光系统的平滑化。利用前端部分相干散射源和严格...
  • 在惯性约束聚变中,辐照的均匀性直接影响到内爆实验的效果,通常结合各种空间、时间平滑技术提高光束的均匀性。光谱色散平滑技术正是一种常用的时间平滑技术,其核心是频谱展宽技术。为提高频谱展宽的效果,将单通的腔式...
  • 图像处理-空间平滑滤波

    千次阅读 2019-08-16 12:51:26
    个人博客:http://www.chenjianqu.com/ 原文链接:... 目录: 1.空间域和频域的概念 2.图像滤波 3.图像卷积 4.常用空间域滤波器 5.线性点运算 6.均值滤波 7.加权均值滤波 ...

    个人博客:http://www.chenjianqu.com/

    原文链接:http://www.chenjianqu.com/show-12.html

    目录:

    1.       空间域和频域的概念

    2.       图像滤波

    3.       图像卷积

    4.       常用空间域滤波器

    5.       线性点运算

    6.       均值滤波

    7.       加权均值滤波

    8.       高斯滤波

    9.       阈值平均滤波

     

    空间域和频域的概念

    空间域与频率域为我们提供了不同的视角。在空间域中,函数自变量(x,y)被视为二维空间中的一个点,数字图像f(x,y)即为一个定义在二维空间中的矩形区域上的离散函数;换一个角度,如果将f(x,y)视为幅值变化的二维信号,则可以通过某些变换手段(如傅里叶变换、离散余弦变换、沃尔什变换和小波变换等)在频域下对图像进行处理了  因为在频率域就是一些特性比较突出,容易处理。比如在空间图像里不好找出噪声的模式,如果变换到频率域,则比较好找出噪声的模式,并能更容易的处理。

    空间域 英文: spatial domain。 释义: 又称图像空间(image space)。由图像像元组成的空间。在图像空间中以长度(距离)为自变量直接对像元值进行处理称为空间域处理。

    频率域。 英文: spatial frequency domain。 释义: 以频率(即波数)为自变量描述图像的特征,可以将一幅图像像元值在空间上的变化分解为具有不同振幅、空间频率和相位的简振函数的线性叠加,图像中各种频率成分的组成和分布称为空间频谱。这种对图像的频率特征进行分解、处理和分析称为频率域处理或波数域处理。

     

    图像滤波

    图像滤波,即在尽量保留图像细节特征的条件下对目标图像的噪声进行抑制,是图像预处理中不可缺少的操作,其处理效果的好坏将直接响到后续图像处理和分析的有效性和可靠性。(滤波就是要去除没用的信息,保留有用的信息,可能是低频,也可能是高频)。

    滤波的目:1. 是抽出对象的特征作为图像识别的特征模式;  2. 为适应图像处理的要求,消除图像数字化时所混入的噪声。

    滤波可分为空间域滤波可频域滤波,前者通过图像卷积运算实现,后者通过傅立叶变换实现。

     

    图像卷积   

        卷积运算:可看作是加权求和的过程,使用到的图像区域中的每个像素分别于卷积核(权矩阵)的每个元素对应相乘,所有乘积之和作为区域中心像素的新值。

    卷积核:卷积时使用到的权用一个矩阵表示,该矩阵是一个权矩阵。不同的卷积核对应不同的滤波器。

    边界问题:当处理图像边界像素时,卷积核与图像使用区域不能匹配,卷积核的中心与边界像素点对应,卷积运算将出现问题。解决办法:1.复制边界,2.忽略边界,3.设定边界值。

    卷积示例:

    3 * 3 的像素区域R与卷积核G的卷积运算

    1.jpg

    R5(中心像素)=R1G1 + R2G2 + R3G3 + R4G4 + R5G5 + R6G6 + R7G7 + R8G8 + R9G9

    代码实现:(对边缘不进行处理)

    //kernel为float类型的卷积核
    Mat ConvolutionOperation(Mat &src, Mat &kernel)
    {
        Mat dst(src.rows, src.cols, src.type(), Scalar(0));
        if (src.channels() == 1)
        {
            int rowsSub = int(kernel.rows / 2);
            int colsSub = int(kernel.cols / 2);
            for (int i = 0; i < src.rows; i++) {
                for (int j = 0; j < src.cols; j++)
                {
                    if (i < rowsSub || i >= src.rows - rowsSub || j < colsSub || j >= src.cols - colsSub) {
                        dst.at<uchar>(i, j) = src.at<uchar>(i, j);
                    }
                    else {
                        float sum = 0;
                        for (int ki = 0; ki < kernel.rows; ki++)
                        {
                            for (int kj = 0; kj < kernel.cols; kj++)
                            {
                                int i_ = i + ki - rowsSub;
                                int j_ = j + kj - colsSub;
                                sum += src.at<uchar>(i_, j_)*kernel.at<float>(ki, kj);
                            }
                        }
                        dst.at<uchar>(i, j) = int(sum);
                    }
                }
            }
     
        }
        return dst;
    }

    常用的空间域滤波器

    根据滤波的效果,空间域滤波可分为空间平滑滤波器和空间锐化滤波器。平滑滤波报告均值滤波、加权均值滤波、阈值平均滤波、中值滤波、高斯滤波等,应用时他们仅是卷积核之间的不同。

    平滑滤波用于模糊处理和降低噪声。模糊处理常用于预处理任务中,如在目标提取之前去除图像中的一些琐碎细节,以及桥接直线或曲线的缝隙。通过线性或非线性平滑滤波也可降低噪声。

    以下分别介绍各个滤波器。

     

    线性点运算

    对图像进行点运算
    Mat LinearPointOperation_Float(Mat& src, double a, double b)
    {
        Mat dst(src.rows, src.cols, src.type(), Scalar(0));
        if (src.channels() == 1)
            for (int i = 0; i < src.rows; i++) {
                for (int j = 0; j < src.cols; j++)
                    dst.at<float>(i, j) = src.at<float>(i, j)*a+b;
        return dst;
    }

    空间域平滑滤波

    1均值滤波

    卷积核:

    2.png

    代码实现:

    Mat MeanFiltering(Mat &src, int n)
    {
        Mat one = Mat::ones(Size(9, 9), CV_32FC1);
        Mat kernel = LinearPointOperation_Float(one, 1 / 81.0, 0);
        Mat dst = ConvolutionOperation(src, kernel);
        return dst;
    }

    运行结果:

    3.png

     

    2.加权均值滤波

    卷积核:

    4.png

    代码实现:

    Mat WeightedMeanFiltering(Mat &src, Mat &kernel)
    {
        Mat dst = ConvolutionOperation(src, kernel);
        return dst;
    }
    Mat WeightedKernel = (Mat_<float>(3, 3) << 1, 2, 1, 2, 4, 2, 1, 2, 1);   
    Mat kernel = LinearPointOperation_Float(WeightedKernel, 1 / 16.0, 0);
    Mat dst = WeightedMeanFiltering(img, kernel);

    运行效果:

    5.png

     

    3. 高斯滤波器

    一维高斯分布:

    6.png

    二维高斯分布:

    7.png

    高斯模板:

    8.png

    9.png

     

    代码实现:

    Mat GaussFiltering(Mat &src, Mat &kernel)
    {
       Mat dst = ConvolutionOperation(src, kernel);
       return dst;
    }
    Mat GaussKernel = (Mat_<float>(5, 5) << 1, 4, 7, 4, 1, 4, 16, 26, 16,4,7,26,41,26,7,4,16,26,16,4,1,4,7,4,1);
    Mat kernel = LinearPointOperation_Float(GaussKernel, 1 / 273.0, 0);
    Mat dst = GaussFiltering(img, kernel);
    imshow("dst", dst);

    运行结果:

    10.png

     

    4. 阈值平均滤波

    顾名思义,阈值平均滤波就是在平均滤波的基础上加上阈值的约束,即当像素点与图像均值的差小于设定的阈值时,输出像素点等于原像素点,否则求其对该像素点进行求模板均值。

    代码实现:

    求图像均值:

    double ImageMean(Mat &src)
    {
        long sum = 0;
        if (src.channels() == 1) {
            for (int i = 0; i < src.rows; i++) 
            {
                uchar *srcRow = src.ptr(i);
                for (int j = 0; j < src.cols; j++)
                    sum += srcRow[j];
            }
        }
        return double(sum*1.0 / (src.cols*src.rows));
    }

    实现阈值均值滤波:

    Mat ThresholdMeanFiltering(Mat &src,int n,int thre)
    {
        //获得卷积核
        Mat kernel = Mat::ones(Size(n, n), CV_32FC1);
        kernel = LinearPointOperation_Float(kernel, 1 / (1.0*n*n), 0);
        //创建输出图像
        Mat dst(src.rows, src.cols, src.type(), Scalar(0));
        if (src.channels() == 1)
        {
            double m = ImageMean(src);//计算出图像均值
            int rowsSub = int(kernel.rows / 2);
            int colsSub = int(kernel.cols / 2);
            for (int i = 0; i < src.rows; i++) {
                for (int j = 0; j < src.cols; j++)
                {
                    if (i < rowsSub || i >= src.rows - rowsSub || j < colsSub || j >= src.cols - colsSub) {
                        dst.at<uchar>(i, j) = src.at<uchar>(i, j);
                    }
                    else {
                        if (abs(src.at<uchar>(i, j)-m) < thre)
                            dst.at<uchar>(i, j) = src.at<uchar>(i, j);
                        else {
                            float sum = 0;
                            for (int ki = 0; ki < kernel.rows; ki++)
                            {
                                for (int kj = 0; kj < kernel.cols; kj++)
                                {
                                    int i_ = i + ki - rowsSub;
                                    int j_ = j + kj - colsSub;
                                    sum += src.at<uchar>(i_, j_)*kernel.at<float>(ki, kj);
                                }
                            }
                            dst.at<uchar>(i, j) = int(sum);
                        }
                    }
                }
            }
     
        }
        return dst;
    }

    当输入n=5,thre=50时,运行结果:

    11.png

     

    参考文献

    [1] CSDN博客:yeler082. 图像处理技术上的空间域和空间频率域.

    https://blog.csdn.net/yeler082/article/details/78374818. 2017-10-28

    [2] CSDN博客:John9ML. 图像处理基本概念——卷积,滤波,平滑.

    https://blog.csdn.net/weixin_38570251/article/details/82054185. 2018-08-25.

    [3] 韩九强,杨磊.数字图像处理.西安交通大学出版社.2018-08

     

     

     

     

     

     

    展开全文
  • 土壤水分空间插值的克里金平滑效应修正方法.pdf
  • 基于空间平滑MUSIC算法的相干信号DOA估计(2)

    千次阅读 多人点赞 2020-10-31 11:06:09
    空间平滑MUSIC算法(2) 继续上一篇博客,继续讲后向空间平滑和前/后向空间平滑MUSIC算法。 基于空间平滑MUSIC算法的相干信号DOA估计(1) 2.3 后向空间平滑算法 后向空间平滑更准确的说是共轭后向空间平滑,它是对后...

    空间平滑MUSIC算法(2)

    继续上一篇博客,继续讲后向空间平滑和前/后向空间平滑MUSIC算法。
    基于空间平滑MUSIC算法的相干信号DOA估计(1)

    2.3 后向空间平滑算法

    后向空间平滑
    后向空间平滑更准确的说是共轭后向空间平滑,它是对后向子阵列地共轭接收数据协方差矩阵进行平滑。定义第一个共轭后向子阵列由 { M , M − 1 , ⋯   , M − p + 1 } \{M, M-1, \cdots, M-p+1\} {M,M1,,Mp+1} 组成, 第二个子阵列由 { M − 1 , M − 2 , ⋯   , M − p } \{M-1, M-2, \cdots, M-p\} {M1,M2,,Mp} 组成,依次组成的子阵列个数为 L = M − p + 1 个 \mathrm{L}=\mathrm{M}-\mathrm{p}+1 个 L=Mp+1
    容易知道,共轭后向空间平滑协方差矩阵 R ~ b \tilde{R}^{b} R~b与前向空间平滑协方差矩阵 R ~ f \tilde{R}^{f} R~f的关系:
    R ~ b = J ∗ ( R ~ f ) ∗ ∗ J \tilde{R}^{b}=J*(\tilde{R}^{f})^**J R~b=J(R~f)J利用后向空间平滑协方差矩阵和MUSIC算法也可以分辨出多个相干信号的方位。可以证明,该方法最大也可以检测M/2个相干的信号。

    MATLAB代码

    % bss.m
    % Code For Music Algorithm Based On Backward Spatial Spectrum
    % The Signals Are All Coherent
    % Author:痒羊羊
    % Date:2020/10/29
    
    clc; clear all; close all;
    %% -------------------------initialization-------------------------
    f = 500;                                        % frequency
    c = 1500;                                       % speed sound
    lambda = c/f;                                   % wavelength
    d = lambda/2;                                   % array element spacing
    M = 20;                                         % number of array elements
    N = 100;                                        % number of snapshot
    K = 6;                                          % number of sources
    L = 10;                                         % number of subarray
    L_N = M-L+1;                                    % number of array elements in each subarray
    coef = [1; exp(1i*pi/6);... 
            exp(1i*pi/3); exp(1i*pi/2);... 
            exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
    doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals
    
    %% generate signal
    dd = (0:M-1)'*d;                                % distance between array elements and reference element
    A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
    S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
    X = A*(coef*S);                                 % received data without noise, M*N
    X = awgn(X,10,'measured');                      % received data with SNR 10dB
    
    %% reconstruct convariance matrix
    %% calculate the covariance matrix of received data and do eigenvalue decomposition
    Rxx = X*X'/N;                                   % origin covariance matrix
    H = fliplr(eye(M));                             % transpose matrix
    Rxxb = H*(conj(Rxx))*H;
    Rf = zeros(L_N, L_N);                           % reconstructed covariance matrix
    for i = 1:L
        Rf = Rf+Rxxb(i:i+L_N-1,i:i+L_N-1);
    end
    Rf = Rf/L;
    [U,V] = eig(Rf);                                % eigenvalue decomposition
    V = diag(V);                                    % vectorize eigenvalue matrix
    [V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
    U = U(:,idx);                                   % reset the eigenvector
    P = sum(V);                                     % power of received data
    P_cum = cumsum(V);                              % cumsum of V
    
    %% define the noise space
    J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
    J = J(1);                                       % number of principal component
    Un = U(:,J+1:end);
    
    %% music for doa; seek the peek
    dd1 = (0:L_N-1)'*d;
    theta = -90:0.1:90;                             % steer theta
    doa_a = exp(-1i*2*pi*dd1*sind(theta)/lambda);   % manifold array for seeking peak
    music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
    music = 10*log10(music/max(music));             % normalize the result and convert it to dB
    
    %% plot
    figure;
    plot(theta, music, 'linewidth', 2);
    title('Music Algorithm For Doa', 'fontsize', 16);
    xlabel('Theta(°)', 'fontsize', 16);
    ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
    grid on;
    

    后向空间平滑可以看到,在6个入射信号均相干的情况下,基于后向空间平滑的MUSIC算法能较好地对其进行DOA估计,且估计精度较高。

    2.3 前/后向空间平滑算法

    把前向和共轭后向空间平滑协方差矩阵定义为前向空间平滑协方差矩阵和共轭后向空间平滑协方差矩阵的平均值,即:
    R ‾ = 1 2 ( R ~ f + R ~ b ) \overline{\mathbf{R}}=\frac{1}{2}\left(\tilde\mathbf{R}^{f}+\tilde\mathbf{R}^{b}\right) R=21(R~f+R~b)那么只要空间平滑的次数大于等于相干信号源的个数,前向和共轭后向空间平滑协方差矩阵一般情况下都是满秩的。使用前/后向空间平滑的方法最多可以检测的相干信号源个数为 2 M / 3 2 \mathrm{M} / 3 2M/3 个。你可能很好奇这个最大相干信号源检测数是怎么得到的?
    假设: 阵列天线的阵元数为 M \mathrm M M个,前/后向空间平滑次数分别为 L \mathrm L L次,则每个子阵的阵元数为 N = M − L + 1 \mathrm N={\mathrm M}-\mathrm{L}+1 N=ML+1 个, 同时可以知道, 可以分辨的最大信号个数为 M − L \mathrm{M}-\mathrm{L} ML个, 也就是子阵的阵元个数减 1 ; 1 ; 1;前后向分别平滑 N \mathrm{N} N次可以解相干信号的个数为 2 L 2 \mathrm{L} 2L个, 最大情况下,两者相等所以 M − L = 2 L \mathrm M- \mathrm L=2 \mathrm L ML=2L,也就是 L = M / 3 ; \mathrm{L}=\mathrm{M} / 3 ; L=M/3; 所以 2 L = 2 M / 3 , 2 \mathrm{L}=2 \mathrm{M} / 3, 2L=2M/3, 所以前/后向空间平滑最大可以解相干的信号个数为 2 M / 3 2 \mathrm{M} / 3 2M/3 个。所以说采用前/后向空间平滑的改进技术可以很大地提高阵列孔径。

    MATLAB代码

    % fbss.m
    % Code For Music Algorithm Based On Forward And Backward Spatial Spectrum
    % The Signals Are All Coherent
    % Author:痒羊羊
    % Date:2020/10/29
    
    clc; clear all; close all;
    %% -------------------------initialization-------------------------
    f = 500;                                        % frequency
    c = 1500;                                       % speed sound
    lambda = c/f;                                   % wavelength
    d = lambda/2;                                   % array element spacing
    M = 20;                                         % number of array elements
    N = 100;                                        % number of snapshot
    K = 6;                                          % number of sources
    L = 10;                                         % number of subarray
    L_N = M-L+1;                                    % number of array elements in each subarray
    coef = [1; exp(1i*pi/6);... 
            exp(1i*pi/3); exp(1i*pi/2);... 
            exp(2i*pi/3); exp(1i*2*pi)];            % coherence coefficient, K*1
    doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals
    
    %% generate signal
    dd = (0:M-1)'*d;                                % distance between array elements and reference element
    A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
    S = sqrt(2)\(randn(1,N)+1i*randn(1,N));         % vector of random signal, 1*N
    X = A*(coef*S);                                 % received data without noise, M*N
    X = awgn(X,10,'measured');                      % received data with SNR 10dB
    
    %% reconstruct convariance matrix
    %% calculate the covariance matrix of received data and do eigenvalue decomposition
    Rxx = X*X'/N;                                   % origin covariance matrix
    H = fliplr(eye(M));                             % transpose matrix
    Rxxb = H*(conj(Rxx))*H;
    Rxxfb = (Rxx+Rxxb)/2;
    Rf = zeros(L_N, L_N);                           % reconstructed covariance matrix
    for i = 1:L
        Rf = Rf+Rxxfb(i:i+L_N-1,i:i+L_N-1);
    end
    Rf = Rf/L;
    [U,V] = eig(Rf);                                % eigenvalue decomposition
    V = diag(V);                                    % vectorize eigenvalue matrix
    [V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
    U = U(:,idx);                                   % reset the eigenvector
    P = sum(V);                                     % power of received data
    P_cum = cumsum(V);                              % cumsum of V
    
    %% define the noise space
    J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
    J = J(1);                                       % number of principal component
    Un = U(:,J+1:end);
    
    %% music for doa; seek the peek
    dd1 = (0:L_N-1)'*d;
    theta = -90:0.1:90;                             % steer theta
    doa_a = exp(-1i*2*pi*dd1*sind(theta)/lambda);   % manifold array for seeking peak
    music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
    music = 10*log10(music/max(music));             % normalize the result and convert it to dB
    
    %% plot
    figure;
    plot(theta, music, 'linewidth', 2);
    title('Music Algorithm For Doa', 'fontsize', 16);
    xlabel('Theta(°)', 'fontsize', 16);
    ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
    grid on;
    

    前/后向空间平滑因为前/后向空间平滑的改进技术很大地提高了阵列孔径,从以上的DOA结果图可以看出分辨率提高了。

    参考论文:《阵列信号处理相关技术研究,陈四根》
    代码Code均为原创。
    欢迎转载,表明出处。

    展开全文
  • 博客(上)为数字图像处理课程理论,博客(下)为对应的实验部分。 教材: 中文教材:数字图像处理_第三版_冈萨雷斯 ...上一篇博客:空间滤波(上) 链接:https://mp.csdn.net/console/editor/html/...

    博客(上)为数字图像处理课程理论,博客(下)为对应的实验部分。

    教材:
    中文教材:数字图像处理_第三版_冈萨雷斯
    实验教材(matlab版):数字图像处理(MATLAB版)冈萨雷斯 
    英文教材:Digital Image Processing_3ed_Gonzalez

    上一篇博客:空间滤波(上)
    正文:

    空间平滑滤波器


    1. 空间平滑滤波器用途

    作用:用于模糊处理和减少噪声。

    1.1 我们为什么要使用这类滤波器呢?

    原因: 典型的随机噪声由灰度级的急剧变化组成;平滑处理会降低图像的“尖锐”变化;
    “负面效应”:图像边缘模糊化(平滑);(图像边缘也由灰度级的急剧变化组成)

    1.2 空间平滑滤波器中”平滑“俩字的含义

    “滤波”是指接受(通过)或者拒绝一定的频率成分。例如,通过低频的滤波器称为低通滤波器,低通滤波器的最终效果是模糊(平滑)一幅图像。
    在上一篇博客中我们已经讲了,空间域滤波是以一种把模版运算运用于图像的空间域增强的技术;依据滤波频率空间域滤波分为 平滑滤波(减弱和去除高频分量)和 锐化滤波(减弱和去除低频分量)。

    2. 俩类典型的空间平滑滤波器

    2.1 均值滤波器(线性)

    均值滤波器:用包含在滤波 掩模邻域内的像素的平均灰度值 去代替每个像素点的值。
    均值滤波的模版就是ones(n, n),模版内全部元素均是1,即他们的权重一模一样。

    其它经常使用的线性滤波还有:

    • 加权滤波:通常中心元素权重较大,且对称向外递减;
    • 高斯滤波:加权滤波的特例,依据高斯分布确定模版系数。
    盒滤波器(左)和加权均值滤波器(右)

    均值滤波器的缺点:会使图像变得模糊,原因是它对所有的点都是同等对待,在将噪声点分摊的同时,将景物的边界点也分摊了;
    为了改善效果,我们可以选择采用加权平均的方式来构造滤波器。
    均值滤波器通式

    模板尺寸对滤波效果的影响:

        模板尺寸越大,图像越模糊,丢失得图像细节越多。

    平滑空域滤波的缺点及问题:

        如果我们图像处理的主要目的是去除噪声,那么平滑滤波器在去除噪音的过程中也会钝化图像的边和尖锐部分。

    2.2 中值滤波器(非线性)

    中值滤波器:先将掩模内欲求的像素及其邻域的像素值排序(升序或降序),确定 出中值 , 将中值赋予该像素点。
    主要功能:使拥有不同灰度的点看起来更接近于它的邻近值。
    主要用途:去除“椒盐”噪声。

    二维中值滤波的窗口形状和尺寸对滤波效果影响较大,不同的图像内容和不同的应用要求,往往采用不同的窗口形状和尺寸,常用的二维中值滤波窗口有线状、方形、圆形、十字形以及圆环形等。窗口尺寸一般先用3×3,再取5×5逐渐增大,直至滤波效果满意为止。就一般经验来讲,对于有缓变的较长轮廓线物体的图像,采用方形或圆形窗口为宜。对于包含有尖顶物体的图像,用十字形窗口,而窗口大小则以不超过图像中最小有效物体的尺寸为宜。如果图像中点、线、尖角细节较多,则不宜采用中值滤波。

    中值滤波算法的特征:

        在去除噪音的过程中也会较好的保留边的锐度和图像细节。

    在图像处理中,尽管中值滤波器是使用的最为广泛的统计排序滤波器,但是这并不意味着它是唯一的。同样,可以在排序之后取最大值来代替相应的像素点的灰度值,对应的滤波器称为最大值滤波器;或者在排序之后取最小的像素值来代替相应的像素点的灰度值,对应的滤波器称为最小值滤波器。

    均值滤波和中值滤波非常基础,均值滤波相当于低通滤波,有将图像模糊化的趋势,对椒盐噪声基本无能为力。中值滤波的优点是可以很好的过滤掉椒盐噪声,缺点是易造成图像的不连续性。

    最大值滤波是用窗口内像元的最大值来代替中心像元的亮度值,可以发现图像中的亮点,并消除图像中的“椒”噪声(亮度值小的噪声)。

    最小值滤波是用窗口内像元的最小值来代替中心像元的亮度值,可以发现图像中的暗点,并消除图像中的“盐”噪声(亮度值大的噪声)。

    均值滤波对高斯噪声表现较好,对椒盐噪声表现较差;

    中值滤波对高斯噪声表现较差,对椒盐噪声表现较好。


    公式在此编译不便,所以均换成了图片形式;码字不易,如若您觉得质量还行,请给个你的肯定就是我的动力,后更请多多关注、指教!谢谢~

                                                                                                                                作于2020.04

    展开全文
  • 实验三 图像空间平滑与锐化(Python实现)
  • 增加视点数量并提供平滑运动视差通常需要大量的空间信息,依靠大量的数据信息最终可以实现模拟真实场景效果的目的。全息立体图可用于显示三维离散图像或一组三维的空间数据,具有良好的观测效果。给出了三种实现平滑...
  • 空间平滑是指直接对源图像数据做空间变换以达到平滑的目的。它是一种邻域运算,即输出图像中任何像素的值是根据输入图像中对应像素周围一定邻域内像素的值重新计算得到的。图像平滑也称为模糊或滤波,是图像处理中...
  • 1.加深对图像增强及边缘检测技术的感性认识,应用MATLAB工具箱自带的处理函数或自己编程完成相关的工作,分析处理结果,巩固所学理论知识。 2.熟练掌握空域滤波中常用的平滑和锐化滤波器,针对不同类型和强度的...
  • IPv4向IPv6平滑过渡技术的研究,姚远,黄晓放,随着Internet的飞速发展,目前基于IPv4的互联网在实际应用中越来越暴露出其不足之处如地址空间的日益耗尽、服务质量、网络安全等问题
  • 【图像处理笔记】平滑空间滤波器

    万次阅读 2015-09-20 11:45:25
    平滑空间滤波器是低频增强的空间滤波技术。它的目的有两类:一是模糊处理,二是降低噪声。 本文介绍的平滑空间滤波器也分为两类,一类是线性滤波器,比如最简单的简单平均法。但是大多数线性滤波器具有低通特性,...
  • Capon空间平滑解相干

    2010-01-16 14:44:49
    普通的Capon波束形成不能解相干程序,利用空间平滑技术可以解相干,本程序就是利用空间平滑原理解相干
  • 我国经济的快速发展为私人汽车提供了巨大的发展空间,对私家车保有量进行科学预测是确定我国公路交通长短期发展规划,解决交通拥堵、能源紧张、环境污染、资源短缺等一系列问题的关键。首先根据我国2000-2011年的...
  • 【Python CUDA版】河北工业大学数字图像处理实验三:空间平滑与锐化
  • 参数空间过大 数据稀疏严重,词同时出现的情况可能没有,组合阶数高时尤其明显。 2. n-gram模型 为了解决第一个问题引入马尔科夫假设(Markov Assumption):一个词的出现仅与它之前的若干个词有关: 然后利用极...
  • 可以通过适当学习的平滑方法来缓解稳健的过度拟合 本文的代码 陈天龙*,张振宇*,刘思佳,常世宇,... 平整崎input的输入空间 先决条件 pytorch 1.5.1 火炬视觉0.6.1 advertorch 0.2.3 用法 标准培训: python -u

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 49,986
精华内容 19,994
关键字:

空间平滑技术