（包含代码）利用机器学习进行分步信号处理：主成分分析（PCA）、独立成分分析（ICA）与非负矩阵分解（NMF）

def pca(img, n_components=9):    """    Perform PCA dimensionality reduction on the input image    :param img: a numpy array of shape (n_samples, dim_images)    :param n_components: number of principal components for projection    :return: image in PCA projection, a numpy array of shape (n_samples, n_components)    """
    # Compute the covariance matrix    cov_mat = np.cov(img.T)
    # Compute the eigenvectors and eigenvalues    eig_val_cov, eig_vec_cov = np.linalg.eig(cov_mat)
    # Make a list of (eigenvalue, eigenvector) tuples    eig_pairs = [        (np.abs(eig_val_cov[i]), eig_vec_cov[:, i]) for i in range(len(eig_val_cov))    ]

  # Select n_components eigenvectors with largest eigenvalues, obtain subspace transform matrix    eig_pairs.sort(key=lambda x: x[0], reverse=True)    eig_pairs = np.array(eig_pairs)    matrix_w = np.hstack(        [eig_pairs[i, 1].reshape(img.shape[1], 1) for i in range(n_components)]    )
    # Return samples in new subspace    return np.dot(img, matrix_w), matrix_w

def inverse_pca(img, components):    """    Obtain the reconstruction after PCA dimensionality reduction    :param img: Reduced image, a numpy array of shape (n_samples, n_components)    :param components: a numpy array of size (original_dimension, n_components)    :return:    """    reconstruct = np.dot(img, components.T).astype(int)    return reconstruct.reshape(-1,28,28)    # Perform PCA transformpca_images, components = pca(flat_images, k)
# Check the dimension of the reduced datasetprint(pca_images.shape)
# Reconstruct from PCA reductionreconstruct = inverse_pca(pca_images, components)

def ica(X, step_size=1, tol=1e-8, max_iter=10000, n_sources=4):    m, n = X.shape
    # Initialize random weights    W = np.random.rand(n_sources, m)        for c in range(n_sources):        # Copy weights associated to the component and normalize        w = W[c, :].copy().reshape(m, 1)        w = w / np.sqrt((w ** 2).sum())

      for i in range(max_iter):            # Dot product of weight and input            v = np.dot(w.T, X)
            # Pass w*s into contrast function g            gv = np.tanh(v*step_size).T
            # Pass w*s into g prime            gdv = (1 - np.square(np.tanh(v))) * step_size
            # Update weights            wNew = (X * gv.T).mean(axis=1) - gdv.mean() * w.squeeze()

 # Decorrelate and normalize weights            wNew = wNew - np.dot(np.dot(wNew, W[:c].T), W[:c])            wNew = wNew / np.sqrt((wNew ** 2).sum())
            # Calculate limit condition            lim = np.abs(np.abs((wNew * w).sum()) - 1)            # Update weights            w = wNew            # Check for convergence            if lim < tol:                break
        W[c, :] = w.T    return W

def nmf(X, tol=1e-6, max_iter=5000, n_components=4):
    # Random initialization of weights    W = np.random.rand(X.shape[0], n_components)    H = np.random.rand(n_components, X.shape[1])
    # A very big number for initial error    oldlim = 1e9
    # A very small number to ensure matrices are strictly positive    eps = 1e-7

for i in range(max_iter):        # Multiplicative update steps        H = H * ((W.T.dot(X) + eps) / (W.T.dot(W).dot(H) + eps))        W = W * ((X.dot(H.T) + eps) / (W.dot(H.dot(H.T)) + eps))
        # Frobenius distance between WH and X        lim = np.linalg.norm(X-W.dot(H), 'fro')
        # Check for convergence        if abs(oldlim - lim) < tol:            break
        oldlim = lim