个人博客：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的卷积运算

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均值滤波

卷积核：

代码实现：

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;
}

运行结果：

 

2.加权均值滤波

卷积核：

代码实现：

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);

运行效果：

 

3. 高斯滤波器

一维高斯分布：

二维高斯分布：

高斯模板：

 

代码实现：

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);

运行结果：

 

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时，运行结果：

 

参考文献

[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

 

 

 

 

 

 

空间平滑MUSIC算法（2）

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

2.3 后向空间平滑算法

后向空间平滑更准确的说是共轭后向空间平滑，它是对后向子阵列地共轭接收数据协方差矩阵进行平滑。定义第一个共轭后向子阵列由 { M , M − 1 , ⋯   , M − p + 1 } \{M, M-1, \cdots, M-p+1\} {M,M−1,⋯,M−p+1} 组成, 第二个子阵列由 { M − 1 , M − 2 , ⋯   , M − p } \{M-1, M-2, \cdots, M-p\} {M−1,M−2,⋯,M−p} 组成，依次组成的子阵列个数为 $\mathrm{L}=\mathrm{M}-\mathrm{p}+1个$ 。
容易知道，共轭后向空间平滑协方差矩阵${\tilde{R}}^b$与前向空间平滑协方差矩阵${\tilde{R}}^f$的关系：
$${\tilde{R}}^b=J\ast ({\tilde{R}}^f{)}^{\ast }\ast 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 前/后向空间平滑算法

把前向和共轭后向空间平滑协方差矩阵定义为前向空间平滑协方差矩阵和共轭后向空间平滑协方差矩阵的平均值，即：
$$\overline{\mathbf{R}}=\frac{1}{2}\left({\tilde{\mathbf{R}}}^f+{\tilde{\mathbf{R}}}^b\right)$$那么只要空间平滑的次数大于等于相干信号源的个数，前向和共轭后向空间平滑协方差矩阵一般情况下都是满秩的。使用前/后向空间平滑的方法最多可以检测的相干信号源个数为 $2\mathrm{M}/3$ 个。你可能很好奇这个最大相干信号源检测数是怎么得到的？
假设: 阵列天线的阵元数为$\mathrm{M}$个，前/后向空间平滑次数分别为$\mathrm{L}$次，则每个子阵的阵元数为 $\mathrm{N}=\mathrm{M}-\mathrm{L}+1$ 个, 同时可以知道, 可以分辨的最大信号个数为$\mathrm{M}-\mathrm{L}$个, 也就是子阵的阵元个数减$1;$前后向分别平滑$\mathrm{N}$次可以解相干信号的个数为 $2\mathrm{L}$个, 最大情况下，两者相等所以 $\mathrm{M}-\mathrm{L}=2\mathrm{L}$,也就是 $\mathrm{L}=\mathrm{M}/3;$ 所以 $2\mathrm{L}=2\mathrm{M}/3,$ 所以前/后向空间平滑最大可以解相干的信号个数为$2\mathrm{M}/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均为原创。
欢迎转载，表明出处。

空间平滑滤波器

1. 空间平滑滤波器用途

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

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

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

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

2.1 均值滤波器（线性）

2.2 中值滤波器（非线性）

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

中值滤波算法的特征：

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

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

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

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

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

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

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

---

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

                                                                                                                            作于2020.04