SVD - 矩阵变换Python

Question

尝试在Python中计算SVD以找到频谱中最重要的元素,并创建一个仅包含最重要部分的矩阵.

在python我有:

u,s,v = linalg.svd(Pxx, full_matrices=True)

这样可以得到3个矩阵; 其中"s"包含与u,v相对应的大小.

为了构造一个包含信号的所有重要部分的新矩阵,我需要捕获"s"中的最高值并将它们与"u"和"v"中的列匹配,得到的矩阵应该给我数据中最重要的部分.

问题是我不知道如何在Python中执行此操作,例如,如何在"s"中找到最高数字并选择"u"和"v"中的列以创建新矩阵？

(我是Python和numpy的新手)所以任何帮助都将非常感激

编辑:

import wave, struct, numpy as np, matplotlib.mlab as mlab, pylab as pl
from scipy import linalg, mat, dot;
def wavToArr(wavefile):
    w = wave.open(wavefile,"rb")
    p = w.getparams()
    s = w.readframes(p[3])
    w.close()
    sd = np.fromstring(s, np.int16)
    return sd,p

def wavToSpec(wavefile,log=False,norm=False):
    wavArr,wavParams = wavToArr(wavefile)
    print wavParams
    return  mlab.specgram(wavArr, NFFT=256,Fs=wavParams[2],detrend=mlab.detrend_mean,window=mlab.window_hanning,noverlap=128,sides='onesided',scale_by_freq=True)

wavArr,wavParams = wavToArr("wavBat1.wav")

Pxx, freqs, bins = wavToSpec("wavBat1.wav")
Pxx += 0.0001

U, s, Vh = linalg.svd(Pxx, full_matrices=True)
assert np.allclose(Pxx, np.dot(U, np.dot(np.diag(s), Vh)))

s[2:] = 0
new_a = np.dot(U, np.dot(np.diag(s), Vh))
print(new_a)

Answer 1

linalg.svds按降序返回.因此,要选择n最高的数字s,您只需要形成

s[:n]

如果将较小的值设置s为零,

s[n:] = 0

然后矩阵乘法将负责"选择"U和V的适当列.

例如,

import numpy as np
LA = np.linalg

a = np.array([[1, 3, 4], [5, 6, 9], [1, 2, 3], [7, 6, 8]])
print(a)
# [[1 3 4]
#  [5 6 9]
#  [1 2 3]
#  [7 6 8]]
U, s, Vh = LA.svd(a, full_matrices=False)
assert np.allclose(a, np.dot(U, np.dot(np.diag(s), Vh)))

s[2:] = 0
new_a = np.dot(U, np.dot(np.diag(s), Vh))
print(new_a)
# [[ 1.02206755  2.77276308  4.14651336]
#  [ 4.9803474   6.20236935  8.86952026]
#  [ 0.99786077  2.02202837  2.98579698]
#  [ 7.01104783  5.88623677  8.07335002]]

鉴于这里的数据,

import numpy as np
import scipy.linalg as SL
import matplotlib.pyplot as plt

Pxx = np.genfromtxt('mBtasJLD.txt')
U, s, Vh = SL.svd(Pxx, full_matrices=False)
assert np.allclose(Pxx, np.dot(U, np.dot(np.diag(s), Vh)))

s[2:] = 0
new_a = np.dot(U, np.dot(np.diag(s), Vh))
print(new_a)
plt.plot(new_a)
plt.show()

产生