<jupyter_start><jupyter_text>AimImplement turbulence fitting in the frequency domain.<jupyter_code>%load_ext autoreload
%autoreload 2
import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import welch
from scipy import stats, optimize
import sys
from tqdm import tqdm
import warnings

sys.path.append("..")
from src import SystemIdentifier
from src import genpsd
from src.controllers.identifier import log_likelihood
from src.utils import rms
ol = np.load("../data/tt_center_noise/tt_center_noise_nsteps_10000_delay_0.01_dt_21_07_2021_12.npy")
# first, try to replicate Meimon 2010's Figure 7
nsteps = 100000
ident = SystemIdentifier(ol[:2000], fs=1)
# kvals = [0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1, 2, 5, 10]
kvals = [0.001, 0.05, 0.2, 0.4, 0.8, 2, 5]
freqs = np.arange(0.01, 10.0, 0.01)

def make_A(k, f):
    a1 = 2 * np.exp(-2*k*np.pi*f / ident.fs) * np.cos(2*np.pi*f/(ident.fs * np.sqrt((1 - k**2)+0j)))
    a2 = -np.exp(-4*k * np.pi * f / ident.fs)
    return a1, a2
plt.figure(figsize=(8,8))
x0 = [1, 0]
fpeak = 0.01
with np.errstate(invalid="ignore"):
    for k in tqdm(kvals):
        a1, a2 = make_A(k, fpeak)
        states = np.zeros((nsteps,))
        states[:2] = x0
        for i in range(2, nsteps):
            states[i] = a1 * states[i-1] + a2 * states[i-2]# + 1e-4 * np.random.randn()

        f, p = genpsd(states, dt=1, remove_dc=False)
        plt.loglog(f / fpeak, p / p[0], label=k)
    plt.xlabel(r"$f / f_{peak}$")
    plt.ylabel("PSD / PSD(f=0)")
    plt.ylim(1e-4, 1e10)
    plt.legend()<jupyter_output>100%|██████████| 7/7 [00:01<00:00,  3.85it/s]<jupyter_text>So I can do that, cool. I'm convinced that a turbulence pattern can be fit well by the same damped-exp modes we've been using for vibrations. Now to do it.(cells redacted)this doesn't worklet's just fit an AR process<jupyter_code>ident = SystemIdentifier(ol)
klqg = ident.make_klqg_ar(ar_len=4)
klqg.recompute()
from scipy.stats import multivariate_normal as mvn
process_dist = mvn(cov=klqg.W, allow_singular=True)
plt.loglog(*genpsd(klqg.sim_process()[:,0], dt=1e-3))
rms(klqg.sim_process(nsteps=1000))
plt.loglog(f, p)
plt.loglog(*genpsd(klqg.sim_process(nsteps=10000)[:,0], dt=1e-3))
rms(klqg.sim_control(nsteps=1000))
plt.loglog(*genpsd(klqg.sim_control(nsteps=10000)[:,0], dt=1e-3))
p
kl = ident.make_klqg_turb(np.abs(np.fft.ifft(np.fft.fftshift(p))))
f, p = genpsd(ol[:,0], dt=1e-3)
plt.plot(np.abs(np.fft.ifft(np.fft.fftshift(p)[:len(p)//2])))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1718377088305489, 'fraction_numerical': 0.07000795544948289, 'mean_word_length': 4.228690228690229, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '612613', 'license': 'CC0-1.0', 'n_tokens_mistral': 1217, 'n_tokens_neox': 1094, 'n_words': 235, 'path': 'notebooks/turb_fit.ipynb', 'repo_name': 'aditya-sengupta/sealrtc-mirror', 'size': '2514'}