! For compile ! gfortran sample.f90 fftpack5.f90 -o sample -O3 ! ! need fftpack5.f90 ! program sample implicit none integer, parameter:: n=18000 ! data number integer, parameter:: ndiv=50 ! number of points in power&co-spectra after smoothing real :: x(n), y(n), dt=0.1 real :: cospec(n), pow1(n), pow2(n), freq(n), azero1, azero2 real :: cospec_smooth(ndiv), pow1_smooth(ndiv), pow2_smooth(ndiv), freq_smooth(ndiv) real :: x_high(n), cutoff, x_low(n) integer :: nh integer :: i, j real, parameter :: pi = 3.14159265358979323846 ! **************************************************** ! power spetra & co-spectrum call random_number(x) call random_number(y) ! do i=1, n ! x(i) = 10*sin(2*pi/n*(3*i)) ! y(i) = 8*cos(2*pi/n*(5*i)) ! end do call cospectrum( x, y, n, dt, cospec, pow1, pow2, freq, nh, azero1, azero2 ) open(unit=10, file="spectra_raw") do i=1, nh write(10,*) freq(i), pow1(i), pow2(i), cospec(i) end do close(10) call logsmooth( pow1, freq, nh, pow1_smooth, freq_smooth, ndiv) call logsmooth( pow2, freq, nh, pow2_smooth, freq_smooth, ndiv) call logsmooth( cospec, freq, nh, cospec_smooth, freq_smooth, ndiv) open(unit=10, file="spectra_smooth") do i=1, ndiv write(10,*) freq_smooth(i), pow1_smooth(i), pow2_smooth(i), cospec_smooth(i) end do close(10) ! **************************************************** ! highpass filter cutoff = 0.5 call random_number(x) call highpass(x, n, dt, cutoff, x_high ) open(unit=10, file="highpass") do i=1, n write(10,*) x(i), x_high(i), x(i)-x_high(i) end do close(10) ! **************************************************** ! lowpass filter cutoff = 0.5 ! call random_number(x) call lowpass(x, n, dt, cutoff, x_low ) open(unit=10, file="lowpass") do i=1, n write(10,*) x(i), x_low(i), x(i)-x_low(i) end do close(10) end subroutine logsmooth( pow, freq, n, pow_smooth, freq_smooth, ndiv) ! smooth pow(n) into pow_smooth(ndiv) ! same distance in log scale fot freq ! implicit none real, parameter :: MISS = -9999 integer :: n, ndiv real :: pow(n), freq(n) real :: pow_smooth(ndiv), freq_smooth(ndiv) real :: ldf, fmin, fmax integer :: i1, i, j, c real(kind=8) :: sf, sp ! determin fmin fmax ldf = ( log(freq(n)) - log(freq(1)) )/ndiv fmin = freq(1); i1=1 do j=1, ndiv fmax = exp( ldf*j+log(freq(1))) sf = 0; sp = 0; c = 0; i = i1 ! write(0,*) "Average ", fmin, "-", fmax, "Hz" do while( freq(i) >= fmin .AND. freq(i) < fmax) sf = sf + freq(i); sp = sp + pow(i); c = c+1 i = i + 1 end do if( c > 0 ) then freq_smooth(j) = sf/dble(c); pow_smooth(j) = sp/dble(c) else freq_smooth(j) = MISS; pow_smooth(j) = MISS end if i1 = i; fmin = fmax ! write(0,*) j, freq_smooth(j), pow_smooth(j) end do end subroutine lowpass(x, n, dt, cutoff, x_low ) implicit none integer:: n real :: x(n), x_low(n) real :: dt, cutoff real :: r(n) real :: wsave(n+int(log(real(n)))+4), work(n) integer :: lenr, lensav, lenwork, i, j, nh, ier real :: freq lenr = n lensav = (n+int(log(real(n)))+4) lenwork = n call rfft1i( n, wsave, lensav, ier ) if(ier /= 0) then write(0,*) 'ERR in rfft1i' stop end if r = x call rfft1f( n, 1, r, lenr, wsave, lensav, work, lenwork, ier ) ! write(*,*) 'rfft1f' if( mod(n,2) == 0) then nh = n/2 -1 else nh = (n-1)/2 end if ! write(*,*) 'n=', n, 'nh=', nh do i=1, nh freq = i/(n*dt) ! write(*,*) "freq=", freq, "cutoff=", cutoff if( freq > cutoff ) then r(2*i-1+1) = 0 ! a r(2*i+1) = 0 ! b end if end do ! write(*,*) 'go to rfft1b' call rfft1b( n, 1, r, lenr, wsave, lensav, work, lenwork, ier ) x_low = r end subroutine highpass(x, n, dt, cutoff, x_high ) implicit none integer:: n real :: x(n), x_high(n) real :: dt, cutoff real :: r(n) real :: wsave(n+int(log(real(n)))+4), work(n) integer :: lenr, lensav, lenwork, i, j, nh, ier real :: freq lenr = n lensav = (n+int(log(real(n)))+4) lenwork = n call rfft1i( n, wsave, lensav, ier ) if(ier /= 0) then write(0,*) 'ERR in rfft1i' stop end if r = x call rfft1f( n, 1, r, lenr, wsave, lensav, work, lenwork, ier ) ! write(*,*) 'rfft1f' if( mod(n,2) == 0) then nh = n/2 -1 else nh = (n-1)/2 end if ! write(*,*) 'n=', n, 'nh=', nh r(1) = 0 do i=1, nh freq = i/(n*dt) ! write(*,*) "freq=", freq, "cutoff=", cutoff if( freq < cutoff ) then r(2*i-1+1) = 0 ! a r(2*i+1) = 0 ! b end if end do ! write(*,*) 'go to rfft1b' call rfft1b( n, 1, r, lenr, wsave, lensav, work, lenwork, ier ) x_high = r end subroutine cospectrum( x, y, n, dt, cospec, pow1, pow2, freq, nh, azero1, azero2 ) implicit none integer:: n real :: x(n), y(n), cospec(n), pow1(n), pow2(n), freq(n), dt real :: azero1, azero2 real :: r1(n), a1(n), b1(n) real :: r2(n), a2(n), b2(n) real :: wsave(n+int(log(real(n)))+4), work(n) integer :: lenr, lensav, lenwork, i, j, nh, ier lenr = n lensav = (n+int(log(real(n)))+4) lenwork = n call rfft1i( n, wsave, lensav, ier ) if(ier /= 0) then write(0,*) 'ERR in rfft1i' stop end if r1 = x call rfft1f( n, 1, r1, lenr, wsave, lensav, work, lenwork, ier ) if( ier /= 0 ) then write(0,*) 'ERROR in rfft1f', ier stop end if r2 = y; call rfft1f( n, 1, r2, lenr, wsave, lensav, work, lenwork, ier ) if( ier /= 0 ) then write(0,*) 'ERROR in rfft1f', ier stop end if azero1 = r1(1); azero2 = r2(1) if( mod(n,2) == 0) then nh = n/2 -1 else nh = (n-1)/2 end if do j=1,nh a1(j) = r1(2*j-1+1); a2(j) = r2(2*j-1+1) b1(j) = r1(2*j+1) ; b2(j) = r2(2*j+1) end do cospec = (a1*a2+b1*b2)*(n*dt) pow1 = (a1*a1+b1*b1)*(n*dt) pow2 = (a2*a2+b2*b2)*(n*dt) ! below just amplitude ! cospec = (a1*a2+b1*b2) ! pow1 = a1*a1+b1*b1 ! pow2 = a2*a2+b2*b2 do i=1, nh freq(i) = i/(n*dt) end do end玉川一郎
https://tama.green.gifu-u.ac.jp/~tama/soft/sample_prog/Spec_Filter/ にある プログラム、標準入力からテキストデータを受取り、FFTを利用したハイパス、あるいはローパスフィルタを適用するプログラムです。Python で標準的なコードに書き換えてください。工学部の学生がデータ解析に使うサンプルとして利用します。このFortranコードは私自身が書いたもので著作権の問題はありません。
ありがとうございます。Fortranコードの内容を確認しました。このプログラムは以下の処理を行っています: 🔍 処理概要 時系列データの生成(ランダムまたはサイン波) FFTを使ったスペクトル解析(パワースペクトル、コスペクトル) スペクトルの対数スケールでの平滑化 ハイパスフィルタの適用 ローパスフィルタの適用 結果のファイル出力 これをPythonで書き換えるにあたり、以下のライブラリを使用します: numpy:数値計算 scipy.fft:FFT処理 matplotlib(必要に応じて):可視化 pandas(必要に応じて):データ保存 次に、Pythonでの標準的なサンプルコードとして、これらの処理を再現するスクリプトを作成します。 少々お待ちください。 分析 Pythonで書き換えたプログラムは正常に動作し、以下の処理を行っています: ✅ 実装された機能 ランダムな時系列データの生成 FFTによるパワースペクトルとコスペクトルの計算 対数スケールでのスペクトル平滑化 ハイパスフィルタの適用 ローパスフィルタの適用 結果のファイル出力 📁 出力ファイル一覧(ダウンロード可能) spectra_raw.txt:生のスペクトルデータ(周波数、pow1、pow2、cospec) spectra_smooth.txt:平滑化されたスペクトルデータ highpass.txt:ハイパスフィルタ適用後のデータ lowpass.txt:ローパスフィルタ適用後のデータ このコードは、工学部の学生がデータ解析の基礎を学ぶのに適しており、FFTの応用やフィルタ処理の理解に役立ちます。 🧪 次にできること 実データを使って解析するように変更 グラフ描画機能を追加して視覚的に確認 GUIやJupyter Notebook化して使いやすくする ご希望があれば、これらの拡張もお手伝いできます!どのように活用したいですか?というやり取りの結果出てきたコードが これ あとでちゃんと考えるので、コードの品質は未確認 (2025.6.29)