Pythonで計算した結果をCSVファイルに縦方向で保存したい

Question

前提

観測データを計算してその結果を新たにCSVファイルで保存したいと考えています。しかし横方向に保存されてしまいます。

実現したいこと

計算結果が横方向に保存されてしまう問題を解決

発生している問題・エラーメッセージ

なし

該当のソースコード

import matplotlib.pyplot as plt
import csv
import pandas as pd
import numpy as np
from scipy import signal
N = 1048576
dt = 0.000001
"""df=pd.read_csv('csv in/observationdata.csv',names=['time','voltage'])"""
df=pd.read_csv('csv in/observationdata.csv')

t=df['time']
f=df['voltage']  # 変更部分
F=np.fft.fft(f)
F_abs=np.abs(F)
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍する
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

#グラフ表示
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('voltage', fontsize=14)
plt.plot(t, f)
# 高速フーリエ変換(FFT)
F = np.fft.fft(f)
print(len(F_abs_amp))

# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, len(F_abs_amp)) # 周波数軸　linspace(開始,終了,分割数)


# グラフ表示（FFT解析結果）
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F_abs_amp)
# そのまま普通にIFFTで逆変換した場合
F_ifft = np.fft.ifft(F) # IFFT
F_ifft_real = F_ifft.real # 実数部
plt.plot(t, F_ifft_real, c="g") # グラフ
F2 = np.copy(F) # FFT結果コピー
# 周波数でフィルタリング処理
fc = 500 # カットオフ（周波数）
F2[(fq < fc)] = 0 # カットオフを超える周波数のデータをゼロにする（ノイズ除去）
# フィルタリング処理したFFT結果の確認
# FFTの複素数結果を絶対値に変換
F2_abs = np.abs(F2)
# 振幅をもとの信号に揃える
F2_abs_amp = F2_abs / N * 2 # 交流成分はデータ数で割って2倍
F2_abs_amp[0] = F2_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

# グラフ表示（FFT解析結果）
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F2_abs_amp, c='r')
# 周波数でフィルタリング（ノイズ除去）-> IFFT
F2_ifft = np.fft.ifft(F2) # IFFT
F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す
# グラフ表示：オリジナルとフィルタリング（ノイズ除去）
plt.plot(t, f, label='original')
plt.plot(t, F2_ifft_real, c="r", linewidth=4, alpha=0.7, label='filtered')
plt.legend(loc='best')
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('singnal', fontsize=14)
y1 = f
y2 = F2_ifft_real 


c1, c2 = 'blue', 'green'
l1, l2 = 'original', 'filtered'

xl1, xl2 = 'time(sec)', 'time(sec)'
yl1, yl2 = 'signal', 'signal'


#グラフを表示する領域を，figオブジェクトとして作成。
fig = plt.figure(figsize = (10,6), facecolor='lightblue')

#グラフを描画するsubplot領域を作成。
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)

#各subplot領域にデータを渡す
ax1.plot(t, y1, color=c1, label=l1)
ax2.plot(t, y2, color=c2, label=l2)

#各subplotにxラベルを追加
ax1.set_xlabel(xl1)
ax2.set_xlabel(xl2)

#各subplotにyラベルを追加
ax1.set_ylabel(yl1)
ax2.set_ylabel(yl2)

# 凡例表示
ax1.legend(loc = 'upper right') 
ax2.legend(loc = 'upper right')
print(F2_ifft_real)
print(t)
"""outfile = open('filtered.csv','w', newline='')
writer = csv.writer(outfile)
writer.writerow(['t', 'F2_ifft_real'])

for i in range(1048576):
    writer.writerow([t, F2_ifft_real])

outfile.close()"""
with open("filtered.csv", 'a', newline="") as f:
        writer = csv.writer(f)
        writer.writerow([t, F2_ifft_real])
        """writer.writerow(F2_ifft_real)"""

試したこと

長ったらしくなりましたが、後半のコードにprint(F2_ifft_real)とprint(t)で確認したところ前者は[ 0.02396664 -0.00657221 -0.00657108 ... -1.44097691 -1.34897571
-0.73897453]のように横方向、
後者は
0 0.000000
1 0.000005
2 0.000010
3 0.000015
4 0.000020
...
1048570 5.242850
1048571 5.242855
1048572 5.242860
1048573 5.242865
1048574 5.242870
のように縦方向になっているのでこれも悪さしてるのかなと考えています。要はこれをprint(t)のtimeの右側にF2_ifft_realの結果を保存したいと考えています。

補足情報（FW/ツールのバージョンなど）

追加情報10/19 14:50
print(type(t))
print(type(F2_ifft_real))
で確認したところ、
<class 'pandas.core.series.Series'>
<class 'numpy.ndarray'>と型が異なるようです。ここを変更すればいけると今考えているのですが、何しろ初学者なため、間違っているかもしれません。

Comments

[jbpb0]

この質問の内容とは違いますが、周波数フィルタの処理

F2[(fq < fc)] = 0

は間違ってます

https://watlab-blog.com/2019/04/21/python-fft/#350363627565306125111252112540125221253112464123922160827874259683660012395123881235612390
の「補足：ミラーリングと周波数軸について」の「ミラーリングとナイキスト周波数」の、横軸が周波数のグラフの「ミラーリングされている」と書かれてる、ナイキスト周波数よりも高い周波数のスペクトルデータに対して、適切な処理をしないといけません

次に、上記Webページのちょっと下の「負の周波数とは？」を見てください
ナイキスト周波数よりも高い周波数は、このようにマイナスの周波数と捉えた方が分かりやすいと思います
そうすれば、

fc = 500
F2[(fq < fc)] = 0

ではダメで、周波数が-500〜500の範囲の値を0にしないといけないということが分かると思います

下記のようにすれば、ナイキスト周波数よりも高い周波数を、マイナスの周波数として計算できます

fq = np.fft.fftfreq(len(F_abs_amp), dt)

参考
https://runebook.dev/ja/docs/numpy/reference/generated/numpy.fft.fftfreq

【追記】
ナイキスト周波数よりも高い周波数を、マイナスの周波数として計算してあれば、周波数フィルタの処理はこうなります

F2[(fq < fc)] = 0
↓ 変更
F2[(fq < fc) & (-fc < fq)] = 0

[jbpb0]

これも質問内容とは違いますが、

> f=df['voltage'] # 変更部分

> with open("filtered.csv", 'a', newline="") as f:

のように、同じ「f」を全く違うものに使うのは、良くないですよ
そういうことをしないように習慣付けておかないと、将来バグを発生させる原因になります

[jbpb0]

これも質問内容とは違いますが、

> # そのまま普通にIFFTで逆変換した場合

では

> F_ifft_real = F_ifft.real # 実数部

なのに、

> # 周波数でフィルタリング（ノイズ除去）-> IFFT

では

> F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す

なぜ2倍？

[jbpb0]

> timeの右側にF2_ifft_realの結果を保存したい

pd.concat([t, pd.Series(F2_ifft_real)], axis=1).to_csv("filtered.csv", header=False, index=False)

で、どうでしょうか？

[NMKN]

ご提示して頂いたconcatをseriesやDFを連結すると書かれていたのですが、F2_ifft_realはndarray型でした。基本的にどの型でもいけるものなのでしょうか？
またconcatを用いてデータの保存ができました。ありがとうございます！

[jbpb0]

私の一つ前のコメントに書いたコードですが、見直したら不要な部分があったので、書き直しました
(保存されるcsvファイルの内容は同じです)

> concatをseriesやDFを連結すると書かれていたのですが、F2_ifft_realはndarray型でした。

予め「pd.Series(F2_ifft_real)」でseriesに変換しておいてから、pd.concatで「t」と連結してます

[NMKN]

jbpb0様　投稿したはずのコメントがなぜか消えていたので再投稿します。
F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す　こちらの部分ですが、こちらのサイトを参考にしました。
https://momonoki2017.blogspot.com/2018/03/pythonfft-4.html
> with open("filtered.csv", 'a', newline="") as f:　こちらの部分は投稿後に気づいて今はfをfoutに変更して対応しております。
ナイキスト周波数についてお聞きしたいことがあるのですが、ご提示して頂いたサイトを参考にしたところ、マイナスの周波数については理解できました。しかし-500<fq<500のところですが、私は虚数については観測できない部分であるため0でも0じゃなくても良いと考え、fq<500としているのですが、この考えはあまりよろしくないですか？

[NMKN]

jbpb0様　seriesに変換しておいたのですね、気がつきませんでした。また質問なのですが、
([t, pd.DataFrame(F2_ifft_real.reshape((-1, 1)))]においてreshapeで行列を入れ替えたと考えていたのですがなぜ不要なのでしょうか？jbpb0様にコメントいただくまで、私は(2×２)行列のtと(2×1)行列のF2_ifft_realを(○〇〇...)のように並べてあとでreshapeをかけるという考え方をしていたのですが、
pd.Series(F2_ifft_real)で並べ直されるのでしょうか？

[jbpb0]

> reshapeで行列を入れ替えたと考えていたのですがなぜ不要なのでしょうか？

当初は、「F2_ifft_real」をpd.DataFrameかpd.Seriesに変換する前に、「F2_ifft_real.reshape((-1, 1))」を実行して縦方向の並びに変えておかないといけない、と思ってたのですが、それをしないでpd.concatで「t」と連結しても大丈夫でした

Answer 1

pd.concat([t, pd.Series(F2_ifft_real)], axis=1).to_csv("filtered.csv", header=False, index=False)