0

MATLAB で作成した次の関数を Python に変換しようとしています。

function output_phase = fix_phasedata180(phase_data_degrees, averaging_length)

x = exp(sqrt(-1)*phase_data_degrees*2/180*pi);
N = averaging_length;
b = 1/sqrt(N)*ones(1,N);
y = fftfilt(b,x);y = fftfilt(b,y(end:-1:1));y = y(end:-1:1); # This is a quick implementation of filtfilt using fftfilt instead of filter
output_phase = (phase_data_degrees-(round(mod(phase_data_degrees/180*pi-unwrap(angle(y))/2,2*pi)*180/pi/180)*180));
temp = mod(output_phase(1),90);
output_phase = output_phase-output_phase(1)+temp;
output_phase = mod(output_phase,360);
s = find(output_phase>= 180);
output_phase(s) = output_phase(s)-360;

そこで、MATLAB で定義されたこの関数を Python に実装しようとしています。

def fix_phasedata180(data_phase, averaging_length):
    x = np.exp(1j*data_phase*2./180.*np.pi)
    N = averaging_length
    b = 1./np.sqrt(N)*np.ones(N)
    y = fftfilt(b,x)          
    y = fftfilt(b,y[::-1])
    y = y[::-1]
    output_phase = data_phase - np.array(map(round,((data_phase/180.*np.pi-np.unwrap(np.angle(y))/2.)%(2.*np.pi))*180./np.pi/180.))*180
    temp = output_phase[0]%90
    output_phase = output_phase-output_phase[0]+temp
    s = output_phase[output_phase >= 180]
    for s in range(len(output_phase)):
        output_phase[s] = output_phase[s]-360
    return output_phase

関数fftfiltは MATLAB の fftfilt のクローンであると考えていましたが、実行すると次のエラーが発生しました

ValueError                                Traceback (most recent call last)
<ipython-input-40-eb6944fd1053> in <module>()
      4 N = averaging_length
      5 b = 1./np.sqrt(N)*np.ones(N)
----> 6 y = fftfilt(b,x)

D:/folder/fftfilt.pyc in fftfilt(b, x, *n)
     66         k = min([i+N_fft,N_x])
     67         yt = ifft(fft(x[i:il],N_fft)*H,N_fft) # Overlap..
---> 68         y[i:k] = y[i:k] + yt[:k-i]            # and add
     69         i += L
     70     return y

ValueError: could not broadcast input array from shape (0,0) into shape (0)

だから、私の質問は: Python で MATLAB fftfilt に相当するものはありますか? 私の関数の目的はoutput_phase、位相信号の速い変動を修正し、n*90 度のシフトを修正することです。 ここに画像の説明を入力

4

4 に答える 4

0

リンク先の関数は、Matlab 関数に相当する Python です。たまたま壊れただけです。

いずれにせよ、MNEには、fftfilt 関数で使用されるオーバーラップおよび追加メソッドの実装もあります。これはライブラリのプライベート関数であり、Matlab 関数とまったく同等に呼び出すことができるかどうかはわかりませんが、おそらく便利です。ソースコードはhttps://github.com/mne-tools/mne-python/blob/master/mne/filter.py#L41にあります。

于 2016-01-11T15:41:42.343 に答える
0

最後に、コードを改善しました。fftfilt（2回適用）をscipy.signal.filtfilt（基本的に同じ）に置き換えます。したがって、Pythonに翻訳された私のコードは次のようになります。

import numpy as np
import scipy.signal as sg

AveragingLengthAmp = 10
AveragingLengthPhase = 10
PhaseFixLength = 60
averaging_length = channel_sampling_freq1*PhaseFixLength

def fix_phasedata180(data_phase, averaging_length):
    data_phase = np.reshape(data_phase,len(data_phase))
    x = np.exp(1j*data_phase*2./180.*np.pi)
    N = float(averaging_length)
    b, a = sg.butter(10, 1./np.sqrt(N))
    y = sg.filtfilt(b, a, x)
    output_phase = data_phase - np.array(map(round,((data_phase/180*np.pi-np.unwrap(np.angle(y))/2)%(2*np.pi))*180/np.pi/180))*180
    temp = output_phase[0]%90
    output_phase = output_phase-output_phase[0]+temp
    s = output_phase[output_phase >= 180]
    for s in range(len(output_phase)):
        output_phase[s] = output_phase[s]-360
    return output_phase

out1 = fix_phasedata180(data_phase, averaging_length)

def fix_phasedata90(data_phase, averaging_length):
    data_phase = np.reshape(data_phase,len(data_phase))
    x = np.exp(1j*data_phase*4./180.*np.pi)
    N = float(averaging_length)
    b, a = sg.butter(10, 1./np.sqrt(N))
    y = sg.filtfilt(b, a, x)
    output_phase = data_phase - np.array(map(round,((data_phase/180*np.pi-np.unwrap(np.angle(y))/4)%(2*np.pi))*180/np.pi/90))*90
    temp = output_phase[0]%90
    output_phase = output_phase-output_phase[0]+temp
    output_phase = output_phase%360
    s = output_phase[output_phase >= 180]
    for s in range(len(output_phase)):
        output_phase[s] = output_phase[s]-360
    return output_phase

offset = 0
data_phase_unwrapped = np.zeros(len(out2))
data_phase_unwrapped[0] = out2[0]
for jj in range(1,len(out2)):
    if out2[jj]-out2[jj-1] > 180:
        offset = offset + 360
    elif out2[jj]-out2[jj-1] < -180:
        offset = offset - 360
    data_phase_unwrapped[jj] = out2[jj] - offset

fix_phasedata180の場合と同様に、ここで180 度のずれを修正しfix_phasedata90ます。はchannel_sampling_freq11/秒です。

結果は次のとおりです。 ここに画像の説明を入力

それはおおむね正しい。scipy.signal.butterとscipy.signal.filtfiltの理解に疑問があるのは私だけです。ご覧のとおり、私は次のように選択します。

b, a = sg.butter(10, 1./np.sqrt(N))

ここで、フィルターの次数 (N) は 10 で、臨界周波数 (Wn) は 1/sqrt(60) です。私の質問は、フィルターの適切な順序を選択するにはどうすればよいですか? N=1 から N=21 まで試してみましたが、21 より大きい結果data_phase_unwrappedはすべて NAN です。padlen私もin に値を与えてみましたfiltfiltが、よくわかりませんでした。

于 2016-01-18T18:44:42.110 に答える