粗大メモ置き場

個人用,たまーに来訪者を意識する雑記メモ

スキャンしたレシートをGoogle Vision APIを使って自前でOCRしてcsv形式に変換する(Python)

はじめに

概要

レシートをScanSnapでスキャンした画像を
Google Vision APIを使ってOCRして
Pythonを使ってzaimに渡せるようなcsv形式に変換してみた

問題提起

家計管理のため、レシートをzaimという家計管理アプリでスキャンして管理しているのですが自分は定期的に撮影するということができず、200枚近くのレシートをためてしまいました。

ossyaritoori.hatenablog.com

これらの写真を撮影して...という作業に嫌気が指したのでScanSnapでまとめてスキャンしてOCRも自前でできたらいいなということで作業をはじめました。

Google Vision APIについて

Pythonで軽く試せるOCRにはいくつか選択肢があります。
前回Tesseraactを使って見た感想として、チューニングなどしないときちんと精度が出ず結構面倒だったというのがあったので
すでにある程度完成しているGoogle Vision APIというものを使ってみました。

ossyaritoori.hatenablog.com

レシートOCRに関しては下記の記事がちょうど該当したのでこちらを流用してOCRしていこうと思います。

qiita.com

実際の処理

本来は下記のようなデータフローを想定しています。

画像 -> OCR結果  -> 辞書型などのデータ -> csv形式 

しかし、Google Vision APIが従量課金制である都合上OCRをかける回数を最小にしたいため、JSONファイルを介して結果を保存して再利用します。

 (最初だけ行う) 
画像 -> OCR結果  -> JSONファイル

(試行錯誤する処理)
JSONファイル -> 辞書型などのデータ -> csv形式 

画像のOCRと保存

別記事に書いたのでそちらを参照してください。

について書いてあります。

qiita.com

OCR結果を使いやすい形式に変換

下記を流用します。

qiita.com

  • 概要:行ごとにOCR結果のテキストをまとめる
    • 入力:OCR結果のオブジェクト
    • 出力:データのリスト
  • 処理
    • テキストと位置に関する記述を抽出
    • テキストのBoundingBoxの左上の縦方向位置(Y座標)によってテキストをクラスタリング

以降ではこの関数を通して作成したlinesというリストを前提とします。

def get_sorted_lines(response,threshold = 5):
    """Boundingboxの左上の位置を参考に行ごとの文章にParseする

    Args:
        response (_type_): VisionのOCR結果のObject
        threshold (int, optional): 同じ列だと判定するしきい値

    Returns:
        line: list of [x,y,text,symbol.boundingbox]
    """
    # 1. テキスト抽出とソート
    document = response.full_text_annotation
    bounds = []
    for page in document.pages:
        for block in page.blocks:
            for paragraph in block.paragraphs:
                for word in paragraph.words:
                    for symbol in word.symbols: #左上のBBOXの情報をx,yに集約
                        x = symbol.bounding_box.vertices[0].x
                        y = symbol.bounding_box.vertices[0].y
                        text = symbol.text
                        bounds.append([x, y, text, symbol.bounding_box])
    bounds.sort(key=lambda x: x[1])
    # 2. 同じ高さのものをまとめる
    old_y = -1
    line = []
    lines = []
    for bound in bounds:
        x = bound[0]
        y = bound[1]
        if old_y == -1:
            old_y = y
        elif old_y-threshold <= y <= old_y+threshold:
            old_y = y
        else:
            old_y = -1
            line.sort(key=lambda x: x[0])
            lines.append(line)
            line = []
        line.append(bound)
    line.sort(key=lambda x: x[0])
    lines.append(line)
    return lines

必要なテキストの抽出

最低限必要なテキストとして重要な順に下記が挙げられます。

  • 日付
  • 店名
  • 金額
  • 品目名(Optional)

これらについてそれぞれOCR結果を用いて抽出をしていきます。 思ったよりも手順が莫大になってしまったので詳細は個別記事に書きました。

日付の抽出

日付は正規表現を使うことで文字列から抽出できると考えて実装しました。

大まかにはそれで良かったのですが、OCRの都合で年や月といった文字が認識されていなかったりしたのでそういう意味では苦労しました。

qiita.com

店名の抽出

店名は基本的に最初の方にデカデカと載っているため下記の手順を考えました。

  • 上から5行分の文字列を抽出
  • その中から最も縦方向に大きなBoundingBoxを持つものを店名とする

上記は正直必ずしも正しくない方針ですが、ある程度店目に検討が付けばよいと思い強行しました。

qiita.com

金額の抽出

金額は基本的に「小計」または「合計」というキーワードと対応している...と思いましたが思ったよりも様々な表記がなされていて、フローがとても複雑になってしまいました。

qiita.com

OCR結果のcsv形式への変換

さて、最後はZaimにOCR結果をCSVでアップする段ですが、下記のように項目ごとに列を指定できるので普通のcsv形式に変換しちゃって良さそうです。

ファイルアップロード時の画面

Pandasでちょちょっと書けば終わりなのと後述の理由もあってやる気をなくしたので省略させてください。

オチ

散々いろいろやって気づいたのですがZaimは有料会員にならないとカテゴリをcsvから読んでくれません。
ということでこれではカテゴリの管理ができません…。アホかと…。

一応、それっぽくまとめたものを下記においておきます。

github.com

ffmpegを使ってDVDの動画をmp4でPCに保存(要ImgBurn)

はじめに

DVDから動画を抽出したい人けどなんかステマみたいな記事ばかり乱立しているし変なソフトを入れさせられたくない…
ということで手元にあるソフトでなんやかんやした結果できましたよという報告です。

なお,ImgBurnだけはDVDへの書き込み・読み込みどちらでも使えるので許してください。

この界隈は下記の記事のようにフラットに比較や説明をする体で自サイトのソフトを勧めるように誘導するサイトがなんかやたら多く感じます。

dvdcreator.wondershare.jp

用意するもの

  • Windows10
    • ImgBurn : ISOファイルへの変換に必要
    • ffmpeg:mp4への変換に必要

ffmpegもImgBurnもインストールが済んでいるものとします。

手順1:DVD->ISOファイルへと変換

DVDをドライブに差し込んだあと,ImgBurnをつけると下記のようなメニューが出ます。

ImgBurnのメニュー画面

中央左の「Create Image From Disk」をクリックするとDVDの中身をISOファイルへと変換してくれます。

手順2:ISO->mp4変換

下記に詳細があります。

www.kkaneko.jp

差し込んですぐ再生される系のDVD(イベントものでもらえるDVDなど)の場合はコマンドは簡単で

ffmpeg -i <your file>.iso   <output>.mp4

とすれば良いです。

以上!

Adobe Premiere系のソフトが起動しなくなった時にやったこと(Windows10)

サポートコミュニティに投稿すべきな気もしますが質問投稿サイトっぽいレイアウトだったのではてなに放流します。

概要

結論: Creative Cloudをアンインストールして入れ直したら全てうまく行った。


症状と対策

  • Adobe Premiere ProとRushがどちらも起動しなくなった。タスクマネージャーのプロセスにも何故か出て来ない。
  • アンインストールして入れ直そうにもアンインストーラーすらも起動しない。
  • Creative Cloudにも同様の問題があったのでアンインストーラーをダウンロードしてアンインストール
  • Creative Cloud側から諸ソフトを起動するとちゃんと起動した

Adobe公式にはCreative Cloudのアンインストールは慎重に的なことが書いてありますが,今回の件ではまぁしょうがないでしょう。

helpx.adobe.com

参考にしたリンク

まずは公式のリンクを見るべきです。ただ,私のケースではプロセスも立ち上がらなかったので役に立たなかったわけですが。

helpx.adobe.com

露骨なアフィリエイトが入ってはいるが下記の手順もいくつか参考にしました。

0begin.net

ROS2の開発用にWindows Terminalを使ってbash/Powershellのタブ生成時に特定のbatを実行させる話(bashrc的な)

まとめ

  • ROS2 みたいにターミナルを起動後に特定のコマンドを実行しておきたい場合の書き方のメモ
  • ショートカットやWindows Terminalの起動時オプションに特定の記述を書くことで起動時に特定のプロセスを走らせられる
  • 下記のような記法が使える

コマンドプロンプトの場合

# hogefuga.batを実行してそのまま開きっぱなしにする
cmd.exe /k "hogefuga.bat"

Powershellの場合

# cmd1とcmd2を実行してそのまま開きっぱなしにする
powershell.exe -NoExit -Command "cmd1.ps1;  cmd2.ps1"

事例:ROS2でターミナル起動毎に実行する内容

今回の作業目的です。 ROS1でのsource devel/setup.bashに相当する作業がROS2でも必要です。下記チュートリアルを参照。

docs.ros.org

call C:\dev\ros2\local_setup.bat
C:\opt\ros\galactic\x64\local_setup.ps1

上記コマンドをタブを開くたびに入力する必要があり面倒です。 (Powershellにはbashrcみたいなのがあるので実は不要ですが。)

他にもいろいろ設定する必要がありそう?な気がしますが,それは別の記事で検証します。
なんなら一つのファイルにまとめてからそれを実行すれば良いだけの気もしますしね。

補足: ROS2のインストール形式

ROS2のバイナリインストール時にはzipからではなく,aka.ms/rosに示されている手順でインストールしているため若干走らせているコマンドが異なる可能性があります。WindowsにWSLを経由せずにROS2を入れている人の母数が結構少ないのか情報が足りてないので分かり次第追記します。



Windows Terminalを用いた解決策

Windows TerminalはWSLやコマンドプロンプトPowershellなどをまとめて扱えるMS公式のターミナルアプリです。 ここに下図のように呼び出したいターミナルソフトとその起動コマンドを設定することができるため,いちいち実行するのが面倒なコマンドを起動時に処理することができます。

Windows Terminalの設定画面

新しいターミナルの登録手順

ここに新しくROS2用のコマンドプロンプト/Powershellのターミナル設定を登録します。

設定 -> 「新しいプロファイルの追加」から 「コマンドライン」 を使いたいターミナルソフトに応じて編集します。

コマンドプロンプトを用いる時の設定

コマンドプロンプトを用いる場合は下記のように登録すると良いです。(バージョンはgalacticの設定です)

cmd.exe /k "C:\opt\ros\galactic\x64\local_setup.bat"

/k の意味はbatファイルを実行後,プロセスをExitしないという意味です。

追記:その後なんやかんやあってfoxyに変更したのでコマンドは下記に変更しました。

cmd.exe /k "C:\opt\ros\foxy\x64\setup.bat && C:\opt\ros\foxy\x64\share\gazebo\setup.bat"

Powershellを用いるときの設定

Powershellでは下記のような記述を書くことによりROS2の設定を呼び出せます。

powershell.exe -NoExit -Command "C:\opt\ros\galactic\x64\local_setup.ps1"

なお,私はVisualStudioの設定も読ませた方がいいかなと思って下記を登録しています。長い...

powershell.exe -NoExit -Command "&{Import-Module """C:\Program Files (x86)\Microsoft Visual Studio\2019\Community\Common7\Tools\Microsoft.VisualStudio.DevShell.dll"""; Enter-VsDevShell 595a88d4 -SkipAutomaticLocation -DevCmdArguments """-arch=x64 -host_arch=x64""";C:\opt\ros\galactic\x64\local_setup.ps1}"

補足:複数のコマンドを実行したい場合

複数のコマンドを実行する必要性が出た場合は下記のようにすると良いです。

参考にしたURLたち・免責

参考になったりならなかったりしたURLをおいておきます。
筆者が現時点でROS2の開発をフルにやっていないのでこれがおすすめ!とは言えないのであしからず。

windows — Linuxの.bashrcファイルに相当するウィンドウはありますか?

Safest way to run BAT file from Powershell script - Stack Overflow

Start-Process -FilePath """C:\opt\ros\galactic\x64\local_setup.ps1""" -Wait -NoNewWindow みたいな記法が書いてありましたが,これはちょっと期待していたのと動作が違ったのでパス。

3次元回転のオイラー角をそのまま補間したら駄目とは言うが実際どれくらい変わるのか?(Python)

はじめに

3次元回転にはオイラー角、回転行列、クォータニオン、ロドリゲスの回転公式と多岐に渡る表現方法があり、時と場合によってそれらを使い分ける必要があります。

中でも任意の2姿勢を補間する際、「オイラー角を始点と終端で補間するのは悪手で、代わりにLERPやSLERPを使うと良い」というのは教科書などに載っていて皆知っていると思いますが、一体どれほど悪くなるのかについて具体的に示した資料が見つけられなかったので自分でやってみました。

下準備

本題に入る前に下記について作成したものを載せておきます。

  • python/matplotlibを用いた回転角の図示
  • scipyを用いた3D回転の形式変換

matplotlibによる回転角の図示

入力した角度が三次元でどのような見た目になるかについて知りたかったので簡単にmatplotlibで表示できるclassを作りました。

回転順序は 1-2-3(x-y-z)の軸の順に回転する(R = RzRyRx)としました。

# 使い方
plane1  = Plane3D()
plane1.rotate_points_123euler([10,20,20]) # degreeでroll pitch yaw角を入れる。
plane1.show_as_mesh()

こんな感じで回転された平面が図示されます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy
from scipy.spatial.transform import Rotation as R

class Plane3D:
    """
    3次元平面の描画用クラス
    """
    def __init__(self, side_num=10):
        
        assert type(side_num) is int
        self.init_attitude(side_num)
        
        
    def show_as_mesh(self):
        # Create the figure
        fig = plt.figure()
        # Add an axes
        ax = fig.add_subplot(111,projection='3d')
        # plot the surface
        ax.plot_surface(self.__X, self.__Y, self.__Z)
        return ax
    
    def getXYZ(self):
        """
        保持している図示用の点群の座標を吐き出す関数
        """
        return self.__X, self.__Y, self.__Z
        
    def rotate_points_123euler(self,rpy):
        """
        degree単位回転
        1-2-3(x-y-z)の軸の順に回転する。すなわち R = RzRyRxを掛け算する。
        
        """
        R3d = self.euler2Robj(rpy)
        self.normal_vector = R3d.dot(self.normal_vector)
        self.calc_Z()
        
    def euler2Robj(self,rpy):
        """
        degree単位回転行列の取得
        1-2-3(x-y-z)の軸の順に回転する R = RzRyRxを返す
        """
        R3d = R.from_euler('zyx',[rpy[2],rpy[1],rpy[0]],degrees=True).as_matrix()
        return R3d

    def set_points_123euler(self,rpy):
        """
        1-2-3(x-y-z)の軸の順のEuler角に姿勢を固定する
        """
        self.normal_vector = np.array([0.,0.,1.0]).reshape(-1,1)
        self.rotate_points_123euler(rpy)
        
        
    def init_attitude(self,side_num):
        
        # ベースとなるGrid
        self.__X,  self.__Y = np.mgrid[-side_num:side_num+1,-side_num:side_num+1]
        # 法線ベクトル(縦)
        self.normal_vector = np.array([0.,0.,1.0]).reshape(-1,1)
        # 高さ計算
        dlen = side_num*2+1
        self.__Z = np.zeros((dlen,dlen))
            
    def calc_Z(self):
        if self.normal_vector[2,0] != 0:
            Z_ = - self.__X * self.normal_vector[0,0] - self.__Y * self.normal_vector[1,0]
            self.__Z = Z_ /self.normal_vector[2,0]
        else:
            print("Plane is Vertical!!")

ScipyのRotationの使い方について

ScipyのRotationにて1-2-3(x-y-z)の軸の順に回転する R = RzRyRxを返すには下記の書き方をする必要があります。この辺はパッケージによって違うので注意です。

個人的にはr-p-yの順に引数を与えたいと思ったのですが関数と仕様ケースが合わないのでLambda関数などを用いてごまかしています。

# RzRyRxの掛け算(Yaw-Pitch-Roll)
R.from_euler('zyx',[y,p,r]).as_matrix()

ScipyのSLERPの使い方について

ScipyにはSLERP(球面補間)のライブラリが用意されています。 一見すると使い方がわかりにくかったのでこちらもメモしておきます。

from scipy.spatial.transform import Slerp

rpy2ypr = lambda x: [x[2],x[1],x[0]]

# slerp objectを生成
key_rots =R.from_euler('zyx', [rpy2ypr(rpy0), rpy2ypr(rpyn)], degrees = True) # 始点と終点
key_times = [0, 1] # 補間点を設定
slerp = Slerp(key_times, key_rots) # オブジェクト作成

# 0-1の間の補間点における角度をオイラー角で取得
slerp(0.4).as_euler('zyx',degrees=True)

オイラー角補間とSLERP(球面補間)の比較

ここではr-p-yが(0,0,0)と(10,20,30)または(20,40,60)の間で補間をしてみた際の違いを図示します。

SLERPで補間した結果をオイラー角に再度変換し直したものをオイラー角の線形補間と比較すると下図のようになります。

(0,0,0)から(10,20,30)におけるSLERPとオイラー角線形補間の差分

(0,0,0)から(20,40,60)におけるSLERPとオイラー角線形補間の差分

2点間の距離?が離れれば離れるほど差が顕著になるのがわかると思います。

このグラフを見て、結構違うなとなるかあんまり差がないなとなるかはアプリケーション次第だと思います。何か指標のようなものを知っている人がいたら教えて下さい。

ちなみに、もっとも差があるであろう補間の中央付近での3次元形状の違いを図示すると下記の様になります。

(0,0,0)から(20,40,60)の中央点におけるSLERPとオイラー角線形補間の差分の図示

うーん、個人的にはあんまり変わらない気が... 角速度が重要なケースならよりSLERPの重要性が際立ちそうなものですが、以外と動かしてみたらバレないのではというのが筆者の所感です。 なにかCriticalなケースがあれば教えて下さい。

雑談

誤差状態カルマンフィルタ構築に向けたクォータニオンキネマティクス、という良さげな和文テキストが最近公開されました。ちょっと遊んで見ようかなと思っています。

www.flight.t.u-tokyo.ac.jp

DeepLで翻訳したPDFから文章を抽出する(Python,Apache Tika)

はじめに

DeepLの無料版会員は月に3度までPDFやWORD文章などを翻訳にかけることが出来ます。

本記事の概要を下記3行にまとめます。

  • DeepLで翻訳したPDFは保護がかかっており,印刷したり文章をカーソル上でコピーすることができない
  • Python環境があればこれを抽出することができる
  • DeepLの使用規約には抵触していないはずだが指摘されたら消します

参考サイト

Apatch Tikaというドキュメント分析・抽出ツールを使いますが関連する手法などについて下記を参考にしました。

gammasoft.jp

今回はやりませんでしたが直接DeepLのAPIを叩くのも良いと思います。

ai-research-collection.com

抽出手順

環境設定

Windows10のAnaconda3のPython3.6の環境で実行しています。(古くてすみません)

抽出に用いるパッケージは下記コマンドでインストールできました。pipと混ぜるのが嫌な人は別で環境を立てたほうが良いと思います。

pip install tika

抽出と後処理

元のURLを参考にtxtへの書き出しをするプログラムを先に記します。

from tika import parser

file_data = parser.from_file("extract-sample.pdf")  #  この部分でPDFを読み込み
text_ = file_data["content"]                        # テキストを抽出
print(text_)             # 確認

text  = text_.replace("\n\n","")        # 二重改行が多かったので置換して除去

with open("out.txt", "w", encoding="utf-8") as f:
    f.write(text)        # 日本語を扱うのでUnicodeを指定して書き込み

parserが呼ばれている行で変換が行われ,初回実行時にはインストール画面のようなものが出現します。その際には下記のようなログが出ますが気にしなくて良さそうです。

2022-04-24 16:34:20,224 [MainThread  ] [INFO ]  Retrieving http://search.maven.org/remotecontent?filepath=org/apache/tika/tika-server/1.24/tika-server-1.24.jar to C:\Users\xxxx\AppData\Local\Temp\tika-server.jar.
2022-04-24 16:35:14,657 [MainThread  ] [INFO ]  Retrieving http://search.maven.org/remotecontent?filepath=org/apache/tika/tika-server/1.24/tika-server-1.24.jar.md5 to C:\Users\xxxx\AppData\Local\Temp\tika-server.jar.md5.
2022-04-24 16:35:15,777 [MainThread  ] [WARNI]  Failed to see startup log message; retrying...

また,出てきた文書は見た目をベースに改行などのスペーシングも文字として認識されており,私の文書のケースだと二重改行が目立ったので事前にReplaceしておきました。これもPythonを使う利点ですね。

余談:ですます変換

DeepLの翻訳結果は基本的にですます調になり,たまにである調が混ざるという調整が必要なものになります。

WORDにはですます調を検知することはできますが,これを一括で変換してくれないので下記のサイトを使うのが今のところ一番楽そうです。

https://kanasys.com/tech/723

時間波形のsin波のゲイン・位相の変位を計算する(Python/Numpy)

概要

システム同定などのシチュエーションで単一のsin波を入力し、出力される波形とのゲインと位相差を計算したい場面があったので作成しました。

波形がSin波に似通っている&周波数がわかっているならば直交性を用いた解法が良さそうです。

下準備:時系列データ作成

時刻の作成ではnumpyのarangeとlinspaceと仲良くなると良いです。

import numpy as np
import matplotlib.pyplot as plt


N = 1000            # 1000 sample data
sample = 0.001 # 1ms
freq = 5              # 5Hz


# 時刻作成
t = np.arange(0,sample*N,sample)
# 位相ズレ付きのSin波を作成する無名関数
x = lambda phi: np.sin(2*np.pi*freq*t+phi )

# 信号作成
y1 = x(0)
y2 = x(1.3) # 1.3 rad/s ずれたsin波

f:id:ossyaritoori:20220227183859p:plain

直交性を用いた解法

ゲインと位相を知りたい信号が A sin(\omega t + \phi)とすると、 cos(\omega t) sin(\omega t)の直交性から下記のような式が導けます。


F_1(t) = cos(\omega t)A sin(\omega t + \phi) \\
=  Asin(\omega t)cos(\omega t)cos(\phi)+Acos(\omega t)cos(\omega t) sin(\phi) \\
=  \frac{A}{2} sin(2 \omega t) cos(\phi) + \frac{A}{2}(1 + cos(2\omega t) ) sin(\phi)

この  F_1(t) を十分長い時間で平均すると、 \frac{A}{2}sin(\phi)の項のみが残ります。

同様に sin(\omega t)ともとの信号の各時間での値をかけ合わせた数値 F_2(t) を十分長い時間で平均すると  \frac{A}{2}cos(\phi)が残るのでこれらを用いてゲインと位相が計算できる、という手法です。

def get_mag_and_phase(y,x):
    """
    y: input signal
    x: lamda function
    """
    cos_ = x(0)
    sin_ =  x(np.pi/2.0)
    N=len(y)
    a = sin_ * y
    b = cos_ * y
    A,B = a.sum()*2/N ,b.sum()*2/N
    
    mag = np.sqrt(A**2+B**2)
    phase = np.arctan2(A,B) 
    return mag, phase

def calc_mag_and_phase(y1,y2,freq,st):
    """
    input: 
            - y*: signal
            - freq_dot_t: frequency[Hz] * samplingtime [s]
    """
    N = len(y1)
    w = 2* np.pi* freq
    time = np.arange(0,N*st,st)
    x = lambda phi: np.cos(w*time+phi)
    mag1, phase1 = get_mag_and_phase(y1,x)
    mag2, phase2 = get_mag_and_phase(y2,x)
    
    
    return mag2/mag1, phase2-phase1

print(calc_mag_and_phase(y1,y2,freq,sample))

プログラム内では入力と出力の2つの信号の他に一度sinとcosの波を作ってそことの差分からゲインと位相差を計算しています。

これの明確な欠点として、下記の2つが挙げられます。

  • 入出力の信号に無限時間平均で0にならないノイズがのっていると誤差が発生する
  • 入出力信号の長さが周期の整数倍でないと誤差が発生する

前者はRANSACのような繰り返し処理による外れ値除去、後者は周期の整数倍になるように信号の時刻をCropする手が挙げられます。

後者に簡単に対応すると下記のような感じでしょうか。

def calc_mag_and_phase2(y1,y2,freq,st):
    """
    input: 
            - y*: signal
            - freq_dot_t: frequency[Hz] * samplingtime [s]
    """
    N = len(y1)
    w = 2* np.pi* freq
    t_max = np.floor(N*st)
    N_ = int(t_max/st)
    time = np.arange(0,t_max,st)
    x = lambda phi: np.cos(w*time+phi)
    mag1, phase1 = get_mag_and_phase(y1[0:N_],x)
    mag2, phase2 = get_mag_and_phase(y2[0:N_],x)
    
    return mag2/mag1, phase2-phase1

別解

下記サイトに様々な手法が載っています。上記の手法とどう違うのか式を追えてないのですが、

blog.goo.ne.jp

余談:Cross-Correlationで位相計算にずれが生じる

確か波形のズレを調べるだけならCross-Correlationから波形のズレが測れるはずだよな、と思い下記のコードを試してみました。

from scipy import signal

corr = signal.correlate(y1,y2,method='auto')
lags= signal.correlation_lags(len(y1), len(y2))
plt.plot(lags,corr)
amax = corr.argmax()
amin = corr.argmin()

print((lags[amax])*sample*2*np.pi*freq)

結果は 1.2566370614359172 となんとも惜しい値。相互相関の計算では時間軸上で-∞から+∞までの積分を想定しているため繰り返し数が足りないのでは、というお話でした。

stackoverflow.com

他にも位相相関も試してみましたが似たような問題で正しい値が出なかったので注意です。

def phase_corr(sig1,sig2):
    N = len(sig1)
    fft_sig1 = np.fft.fft(sig1)
    fft_sig2 = np.fft.fft(sig2)
    fft_sig2_conj = np.conj(fft_sig2)
    
    R = fft_sig1*fft_sig2_conj
    R/=np.absolute(R)
    r = np.fft.fftshift(np.fft.ifft(R).real)
    ar = np.unravel_index(r.argmax(), r.shape)
    
    plt.plot(r)
    return ar[0]-N/2, r[ar]