生波形同士の位相差を出力したいです．

Question

Pythonで，生波形同士の位相差を出力したいです．

現在，20~30個程度の生波形の入ったcsvファイル（図1を参照）をFFTし，一つのExcelにまとめる処理を行っています．

csvファイルのc列目に生波形1，d列目に生波形2、e列目に生波形3、f列目に生波形4が出力されています。

このCSVファイルをfft解析し、０行目に計測時間を､一行目に生波形1，二行目に生波形2，三行目に生波形3，四行目に生波形6，五行目にサンプリング周波数，六行目にfft解析をした生波形1，七行目にfft解析をした生波形2，八行目にfft解析をした生波形3，九行目にfft解析をした生波形4としています（図２参照）．

今回行いたいこととして，十行目に生波形1と3の位相差 degを，十一行目に生波形1と4の位相差を出力したいです．

つまり、以下のホームページに記載してあるようなことをPython上で行いたいです。
https://blog.goo.ne.jp/je-suis-sumichan/e/0e562abe1dacdfd5ad65763d4658c685
このホームページでいう、信号1Xnを図1の生波形に、信号2Ynを生波形3や4として、処理したいということです。

ホームページ上のstep1.2.3と書かれた列は出力してもしなくてもかまいません。求められた位相差 degという列を出力できれば嬉しいです。

様々なサイトを調べてみたのですが，グラフを出力するものや，コードが膨大なものが多く，既存の下記に示すコードにどのように埋め込めばよいのかわからなくなってしまいました．

下記に図1と図2，既存のコードを記します．
※図2のcsvファイル名と図1のシート名は対応していません．
図2は，入力フォルダ内に1.2.3.4.5.6というcsvが入っていた際にPython を動かしたものとなります．

ご協力お願いいたします．

質問の欄をご覧ください。

図1

図2

import cv2
import numpy as np
import matplotlib.pyplot as plt
import csv
import os
import glob
import pandas as pd


def my_fft(infile):
    ##### 実験条件、分析対象の設定 ######
    t = np.arange(0,8.192,0.0005) #サンプリング時間,周波数の指定(開始時間(=0),終了時間,サンプリング周期)

    file = infile #分析対象のファイル指定

    ##### データの読み込みy #####
    for i in range(4):
        
        def loadvtp(file,i):  #ファイルの読み込みに使用する関数の定義
            vtp = pd.read_csv(file, dtype='float', skiprows=57, usecols=[i+2], nrows=16384, engine='python', header=None)   #解析対象となるデータの読み込み
            vtp = vtp - np.mean(vtp)   #平均を引く
            vtp  = np.squeeze(vtp)
            return vtp

        ##### 主流(x)成分出力の計算 #####

        vx = loadvtp(file,i)

        ##### FFT #####

        freq = np.fft.fftfreq(16384, 1/2000) #周波数軸の生成(サンプリング点数, サンプリング周期)

        f = np.fft.fft(vx) #FFT
        f [0]=f[0]/2 #定常成分を1/2にする(スケーリングに必要)
        amp = np.abs(f)/(16384/2) #FFT結果を振幅に合わせてスケーリング

        freq = freq[0:16384//2] #周波数軸の対称成分を除く
        amp = amp[0:16384//2] #振幅の対称成分を除く

        if i==0:
            yav=vx
            xav=""
            vav=""
            wav=""
            y=amp
            x=""
            v=""
            w=""
        elif i==1:
            yav=yav
            xav=vx
            vav=""
            wav=""
            y=y
            x=amp
            v=""
            w=""
        elif i==2:
            yav=yav
            xav=xav
            vav=vx
            wav=""
            y=y
            x=x
            v=amp
            w=""
        elif i==3:
            yav=yav
            xav=xav
            vav=vav
            wav=vx
            y=y
            x=x
            v=v
            w=amp

    df = [t, yav, xav, vav, wav, freq, y, x, v, w]
    df=pd.DataFrame(df).T

    return df


infile = glob.glob("入力フォルダ名/*.csv")
infile.sort()
print('入力csvファイル数', len(infile))

with pd.ExcelWriter('出力Excel名.xlsx', mode='w') as writer:
    for in_csv_file in infile:
        in_file_name = os.path.basename(in_csv_file)

        fft_df = my_fft(in_csv_file)

        fft_df.to_excel(writer, sheet_name=in_file_name)

Comments

[y_waiwai]

まずはその位相差、というものを定義しよう。
それはどういうものでしょうか。また、どういう結果を想定してるんでしょうか

[nonnon615]

一つ一つのcsvが時系列のようなものになっています。そこで、得られた位相差のアベレージをとり、その時系列での位相差とするつもりです。
そして、横軸に時系列、縦軸に位相差の関係のグラフをプロットすることが目標です。

[nonnon615]

生波形1と、3の位相差がおよそ0となるか180度ずれている、といったような結果がでると予行しており、様々な考察ができるのではないかと考え、このようなプログラムの作成に至りました。

[nonnon615]

定義として、例えば図2を参考にすると、セルc2とe2の位相差をセルL2に、c3とe3の位相差をL3に、といったようにできないか、と考えています。

[nonnon615]

生波形1の列を0行目の時間の列と組み合わせると、横軸が時間、縦軸が生波形1のサインカーブのようなものができます。同じく、生波形2.3.4でも同様にサインカーブがつくれます。それらの位相差が知りたいというわけです。

0秒での位相差、0.0005秒での位相差、と求めれないか、と考えています。

[y_waiwai]

まずはその波形というのをグラフにしてみてはいかがでしょう。
そのうえで、その各波形のグラフ間の、どれを位相差とするのか示しましょう

＃位相差の説明に全くなってませんがな

[nonnon615]

グラフにして、どれを位相差とするかとのことですが、もちろんそれについても考えました。

しかし、グラフにするとゆらぎのあるサイン波のようなものができます。そのため、0.0005秒ごとにどれだけの位相のズレが生じているか求めれないかと考えています

[nonnon615]

説明がわかりにくかったも思いますので、参考としたホームページ等を説明に加えましたので、ご覧ください。

Answer 1

PRIMO信号処理研究所の　Je-suis-Sumichan です。

小生の記事を参照くださりありがとうございます。
実波形から位相差を求めたい、というニーズは多くあるのですが、現在のところ、工学的な方法はかぎられておりそれぞれ一長一短があります。

　私が考案した方法は条件付きなのですが、１）２つの信号の周波数は単一成分、２）振幅は任意の定数　３）それぞれ周波数は一致していなくてもいい　

　特に、単一成分で「波形がずれている＝位相差あり、周波数は同じ」であれば、非常に正確な位相差が計算できます。　ここ一年ほどで、Blogを引用してくださるかたが増え、そろそろ　Python版のサンプルを作ろうと思っています。　計算に際し限られたサンプル数しか必要しません。つまり「限られた有限時間」で計算しており、高い時間分解能を盛っています。しかし限られたサンプル数からの計算、というのは冗長なデータがほとんどなく、ノイズの影響をもろに受けてしまうという、まあ当たり前なことですが、そんな特性をもっています、

　Excel での計算方法をYouTubeに上げてあります。
https://www.youtube.com/watch?v=ZG7wL3z1rHg&t=146s

　周波数が既知で、振幅は任意、位相差を求めたい場合は、下記を参照ください
https://blog.goo.ne.jp/je-suis-sumichan/e/fff3dca22ca958cc1a1a14409ebbfb92
ここまでくると難解なのですが、位相差の計算は「ベクトルの内積・外積」をつかうわかりやすい方法を使っています、
しかしベクトルでやるには信号を直交化しないといけませが、ここで、「正確な瞬時周波数」が必要なのです。この瞬時周波数を求める方法が、重要なところなのです。

　３月中に、上記YouTube動画と同等な計算ができる　python 版を公開したいと思います。

ご参考　
（本方式にて）電力系統周波数を観測した例　https://www.youtube.com/watch?v=ejQsBP8EW80
ほぼ１サイクルごと、周波数をmHz オーダーでリアルタイムの計算し、周波数の急変を検出しました。