压缩感知正交匹跟踪算法（OMP）代码实现之一维连续信号的求解

前天，列举对单一向量OMP算法的求解，这个向量是离散的。

但是我们知道自然界中的信号是模拟（连续）信号。学过一些信号系统的同学应该知道，自然界的信号（光信号，声信号）是连续的。

比如下图的信号

假设要保存上图的信号，我们需要等间隔的采样256次（只是假设，上图可能不需要那么多）。采样256次是奈奎斯特定理规定的最少采样次数，如果低于这个次数，上图就没办法恢复出来了。但是压缩感知却告诉我们，我们可以通过某种方法（测量矩阵）采样，使得采样出来的数据只有64个。这样保存的数据量便从256个下降到64个。这就是压缩感知的迷人之处！！！

这篇代码是在香港大学沙威发布在网上的MATLAB代码的复现。网上我找不到关于OMP算法python3.6的代码，所以试着去复现沙威的代码，一来进一步理解OMP算法，二来可以将这个代码公布出来，方便新手学习。

首先先利用三角函数生成一个信号。如下图所示。

经过一系列计算我们可以近乎完美的恢复出上图信号。恢复出来的信号如下图所示。
好了，直接上代码。

'''
一维连续向量压缩感知的实现（正交匹配跟踪算法）
测量数M>=K*log(N/K),K是稀疏度,N信号长度,可以近乎完全重构
本篇代码是利用香港大学电子工程系的沙威MATLAB代码进行python3.6的复现
编程时间：2019.7.20
编程人：kyda
'''

import matplotlib.pyplot as plt
import numpy as np
import numpy.linalg as lg
from math import sqrt


def get_matri_new_plus(A):
    '''
    求A矩阵的伪逆
    :param A: 输入矩阵
    :return A_new: 求解的伪逆矩阵
    '''
    A_new = A.T.dot(A)
    if A_new.ndim < 2:
        A_new =  1 / A_new
        A_new = A_new * A.T
    else:
        A_new = lg.inv(A_new)
        A_new = A_new.dot(A.T)
    return A_new

def get_A_ba(A):
    '''
    求解A_ba，也就是矩阵列向量单位化
    :param A: 输入矩阵
    :return A_new: 返回列向量单位后的矩阵
    '''
    A_new = A.T
    for row in range(A_new.shape[0]):
        norm2 = lg.norm(A_new[row], ord=2)           # 求向量的2范数
        A_new[row] /= norm2                          # 单位化
    return A_new


def find_pos_VecInMatrix(A, v):
    '''
    求解某个向量在矩阵中的相应行位置
    :param A: 输入矩阵
    :param v: 输入向量
    :return pos: 位置
    '''
    A_rownum = A.shape[0]    # 矩阵行的数量
    for pos in range(A_rownum):
        if (v==A[pos]).all():           # 判断向量v和矩阵中的某一行是否相等
            return pos

'''
1. 时域测试信号生成
'''

K = 7                                           # 稀疏度
N = 256                                         # 信号长度
M = 64                                          # 测量量
f1 = 50                                         # 信号评率1
f2 = 100                                        # 信号频率2
f3 = 200                                        # 信号频率3
f4 = 400                                        # 信号频率4
fs = 800                                        # 采样频率
ts = 1 / fs                                     # 采样间隔
Ts = np.arange(0, N)                            # 采样序列
op = 2 * np.pi * Ts * ts
x = 0.3*np.cos(op*f1) + 0.6*np.cos(op*f2) + 0.1*np.cos(op*f3) + 0.9*np.cos(op*f4)    # 完整信号

'''
2.时域信号压缩传感
'''
Phi = np.random.randn(M, N)                     # 测量矩阵(高斯分布白噪声)
s = Phi.dot(x)                                  # 获得线性测量


'''
3. 正交匹配跟踪算法重构信号(本质上是L_1范数最优化问题)
'''
m = 2 * K                                       # 算法迭代次数(m>=K)
Psi = np.fft.fft(np.eye(N, N)) / sqrt(N)        # 傅里叶正变换矩阵
T = Phi.dot(Psi.T)                              # 恢复矩阵(测量矩阵*正交反变换矩阵)

hat_y = np.zeros(N)                             # 待重构的谱域(变换域)向量
Aug_t = None                                    # 增量矩阵(初始值为空矩阵)
r_n = s                                         # 残差值

product = [0]*N
T_tran = T.T
A_ba = T_tran.copy()                            # 初始化A矩阵
A_ba1 = A_ba                                    # 复制
pos_selected_v = []
for times in range(m):
    w = A_ba.dot(r_n)                           # 求各列向量和残差的投影系数
    pos = np.argmax(abs(w))                     # 最大投影系数所对应的位置
    b_k = A_ba[pos]
    A_ba = np.delete(A_ba, pos, axis=0)         # 选中的列删除
    pos1 = find_pos_VecInMatrix(A_ba1, b_k)
    pos_selected_v.append(pos1)                 # 纪录最大投影系数的位置
    if Aug_t is None:                          # 将最大贡献向量（未单位化的）加入Aug_t矩阵
        Aug_t = T_tran[pos1]
    else:
        Aug_t = np.c_[Aug_t, T_tran[pos1]]
    A_new_plus = get_matri_new_plus(Aug_t)
    aug_y = A_new_plus.dot(s)                   # 解最小二乘解
    r_n = s - Aug_t.dot(aug_y)                  # 求残差


for item in range(m):                           # 重构的谱域向量
    pos = pos_selected_v[item]
    x_r = aug_y[item]
    try:
        hat_y[pos] = x_r
    except np.ComplexWarning:
        pass

hat_x = np.real(Psi.T.dot(hat_y))               # 做逆傅里叶变换重构得到时域信号


plt.plot(hat_x, marker='.', color='r')           # 重建信号
#plt.plot(x)                                     # 原始信号
plt.show()

有问题或者想交流可以添加我的微信公众号