信号重构随笔

import matplotlib.pyplot as plt
import numpy as np
import random
from scipy.linalg import orth

def soft(a,b):
    return np.multiply(np.sign(a), np.maximum(np.abs(a)-b,0))
def ist_basic(y,Phi,lamb,epsilon= 1e-6,loopmax=10000):
    x_in=np.zeros([256,1])
    loop = 0
    while(1):    
        x_in = soft(x_in + np.dot(Phi.T,(y-np.dot(Phi,x_in))),lamb)
        loop = loop + 1
        if(loop>loopmax):
            break
    return x_in     

M = 64    
N = 256      
K = 10   
x =np.zeros([1,N]) 
col =random.sample(range(0,N),K)
Index_K=5*np.random.randn(1,K)
for i in range(len(col)):
    x[0,col[i]]=Index_K[0,i]

Phi = np.random.randn(M,N)
Phi = orth(Phi.T).T
print(Phi.shape)
sigma = 0.005
e = sigma*np.random.randn(M,1)   
y = np.dot(Phi,x.T) + e
lamb = 0.001*np.max(np.abs(np.dot(Phi.T,y)))
x_r = ist_basic(y,Phi,lamb)
plt.figure()
plt.plot(x_r,"-^")
plt.plot(x.T,":o")#原始信号
plt.show()

软阈值迭代算法重构信号的简单例子