
import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
amplitude = 1.5  # 正弦波振幅 [V]
frequency = 30.273e3  # 周波数 [Hz]
sampling_rate = 250e3  # サンプリング周波数 [Hz]
num_points = 256  # サンプリング点数
bit_depth = 8  # 変換bit数

# 入力ノイズの設定
noise_amplitude = 0.1  # ノイズの振幅

# 時間軸
t = np.arange(num_points) / sampling_rate

# 理想的な正弦波信号の生成
pure_signal = amplitude * np.sin(2 * np.pi * frequency * t)

# ノイズの生成
noise = noise_amplitude * np.random.randn(num_points)

# ノイズを含む入力信号の生成
input_signal = pure_signal + noise

# ADC 量子化 (8ビット)
adc_max = 2**(bit_depth - 1) - 1
adc_min = -2**(bit_depth - 1)
quantized_signal = np.clip(np.round(input_signal * adc_max / amplitude), adc_min, adc_max)

# 入力信号のFFT計算
input_fft_values = np.fft.fft(input_signal)
input_fft_magnitude = np.abs(input_fft_values) / num_points
input_fft_magnitude = input_fft_magnitude[:num_points // 2] * 2  # 正の周波数成分のみ

# 量子化信号のFFT計算
fft_values = np.fft.fft(quantized_signal)
fft_magnitude = np.abs(fft_values) / num_points
fft_magnitude = fft_magnitude[:num_points // 2] * 2  # 正の周波数成分のみ

# パワースペクトル計算
input_power_spectrum = (input_fft_magnitude ** 2)  # 入力信号のパワースペクトル
quantized_power_spectrum = (fft_magnitude ** 2)    # 量子化信号のパワースペクトル

# SNDR計算
signal_bin = int(np.round(frequency / (sampling_rate / num_points)))
signal_power = np.sum(fft_magnitude[signal_bin - 1:signal_bin + 2]**2)
noise_power = np.sum(fft_magnitude**2) - signal_power
sndr = 10 * np.log10(signal_power / noise_power)

# ENOB計算
enob = (sndr - 1.76) / 6.02

# SFDR計算
spurious_power = np.max(np.delete(fft_magnitude, signal_bin))
sfdr = 10 * np.log10(signal_power / spurious_power)

# 結果表示
print(f"SNDR: {sndr:.2f} dB")
print(f"ENOB: {enob:.2f} bits")
print(f"SFDR: {sfdr:.2f} dBc")

# FFTプロット設定
freqs = np.fft.fftfreq(num_points, d=1/sampling_rate)

# 波形プロット設定
plt.figure(figsize=(12, 14))

# 時間領域の入力信号（ノイズあり）のプロット
plt.subplot(4, 1, 1)
plt.plot(t, input_signal, label='Input Signal (with noise)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude [V]')
plt.title('Input Signal (with Noise) in Time Domain')
plt.grid(True)
plt.legend()

# 時間領域の量子化信号のプロット
plt.subplot(4, 1, 2)
plt.plot(t, quantized_signal, label='Quantized Signal', color='orange')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude [V]')
plt.title('Quantized Signal in Time Domain')
plt.grid(True)
plt.legend()

# 周波数領域の振幅スペクトルプロット
plt.subplot(4, 1, 3)
plt.plot(freqs[:num_points // 2], 20 * np.log10(input_fft_magnitude), label='Input Signal FFT (with noise)')
plt.plot(freqs[:num_points // 2], 20 * np.log10(fft_magnitude), label='Quantized Signal FFT', color='orange')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude [dB]')
plt.title('FFT of Input and Quantized Signal')
plt.grid(True)
plt.legend()

# パワースペクトルプロット
plt.subplot(4, 1, 4)
plt.plot(freqs[:num_points // 2], 10 * np.log10(input_power_spectrum), label='Input Signal Power Spectrum (with noise)')
plt.plot(freqs[:num_points // 2], 10 * np.log10(quantized_power_spectrum), label='Quantized Signal Power Spectrum', color='orange')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Power [dB]')
plt.title('Power Spectrum of Input and Quantized Signal')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
amplitude = 1.0  # 振幅 [V]（-1～1）
frequency = 1e3  # 周波数 [Hz]
sampling_rate = 20e3  # サンプリング周波数 [Hz]
num_points = 256  # サンプリング点数
bit_depths = [8, 10]  # AD変換ビット数のリスト

# 時間軸
t = np.arange(num_points) / sampling_rate

# Sin波（Ramp波）の生成
input_signal = amplitude * np.sin(2 * np.pi * frequency * t)

# AD変換関数
def adc_quantization(signal, bit_depth, amplitude):
    adc_max = 2**(bit_depth - 1) - 1  # 最大値
    adc_min = -2**(bit_depth - 1)     # 最小値
    # 信号をAD変換し、指定のビット数に量子化
    quantized_signal = np.clip(np.round(signal * adc_max / amplitude), adc_min, adc_max)
    return quantized_signal

# 入力信号と量子化信号のプロット用の設定
fig, ax = plt.subplots(3, 1, figsize=(12, 12))

# 入力信号のプロット
ax[0].plot(t, input_signal, label='Input Signal (Sin wave)')
ax[0].set_xlabel('Time [s]')
ax[0].set_ylabel('Amplitude [V]')
ax[0].set_title('Input Signal in Time Domain')
ax[0].legend()
ax[0].grid(True)

# ビット数ごとにAD変換、FFT、パワースペクトル、SNRを計算
for bit_depth in bit_depths:
    # AD変換（量子化）
    quantized_signal = adc_quantization(input_signal, bit_depth, amplitude)

    # 量子化信号のプロット
    ax[1].plot(t, quantized_signal, label=f'Quantized Signal ({bit_depth}-bit)')
    ax[1].set_xlabel('Time [s]')
    ax[1].set_ylabel('Amplitude [Digital Code]')
    ax[1].set_title('Quantized Signal in Time Domain')
    ax[1].legend()
    ax[1].grid(True)
    
    # FFT計算
    fft_values = np.fft.fft(quantized_signal)
    fft_magnitude = np.abs(fft_values) / num_points
    fft_magnitude = fft_magnitude[:num_points // 2] * 2  # 正の周波数成分のみ

    # パワースペクトル計算
    power_spectrum = fft_magnitude ** 2

    # SN比（SNR）の計算
    signal_bin = int(np.round(frequency / (sampling_rate / num_points)))
    signal_power = np.sum(power_spectrum[signal_bin - 1:signal_bin + 2])  # 信号成分のパワー
    noise_power = np.sum(power_spectrum) - signal_power                   # ノイズ成分のパワー
    snr = 10 * np.log10(signal_power / noise_power)                       # SNRをdBで計算

    # FFTおよびパワースペクトルのプロット
    freqs = np.fft.fftfreq(num_points, d=1/sampling_rate)[:num_points // 2]
    ax[2].plot(freqs, 10 * np.log10(power_spectrum), label=f'Power Spectrum ({bit_depth}-bit) SNR: {snr:.2f} dB')
    ax[2].set_xlabel('Frequency [Hz]')
    ax[2].set_ylabel('Power [dB]')
    ax[2].set_title('Power Spectrum of Quantized Signal')
    ax[2].legend()
    ax[2].grid(True)

# プロットのレイアウト調整と表示
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# サンプリング周波数と信号長の設定
Fs = 1000          # サンプリング周波数 [Hz]
T = 1 / Fs         # サンプリング周期 [秒]
L = 1500           # 信号の長さ
t = np.arange(0, L) * T  # 時間ベクトル

# 信号生成 (振幅 0.8 の DC オフセット、振幅 0.7 の 50 Hz 正弦波、振幅 1 の 120 Hz 正弦波)
S = 0.8 + 0.7 * np.sin(2 * np.pi * 50 * t) + np.sin(2 * np.pi * 120 * t)

# 平均値 0、分散 4 のランダム ノイズを加えた信号
X = S + 2 * np.random.randn(len(t))

# ノイズを含む信号を時間領域でプロット
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(1000 * t, X, linewidth=1)
plt.title("Signal Corrupted with Zero-Mean Random Noise")
plt.xlabel("Time [ms]")
plt.ylabel("X(t)")
plt.grid()

# フーリエ変換の実行
Y = np.fft.fft(X)
P2 = np.abs(Y / L)  # 両側振幅スペクトル
P1 = P2[:L//2 + 1]  # 片側振幅スペクトル
P1[1:-1] = 2 * P1[1:-1]  # 片側振幅スペクトルに変換

# 周波数軸の生成
f = Fs * np.arange(0, (L/2) + 1) / L

# フーリエ変換結果の片側振幅スペクトルをプロット
plt.subplot(3, 1, 2)
plt.plot(f, P1, linewidth=2)
plt.title("Single-Sided Amplitude Spectrum of X(t)")
plt.xlabel("Frequency [Hz]")
plt.ylabel("|P1(f)|")
plt.grid()

# ノイズなしの信号のフーリエ変換
Y_clean = np.fft.fft(S)
P2_clean = np.abs(Y_clean / L)  # 両側振幅スペクトル
P1_clean = P2_clean[:L//2 + 1]  # 片側振幅スペクトル
P1_clean[1:-1] = 2 * P1_clean[1:-1]  # 片側振幅スペクトルに変換

# ノイズなしの信号の片側振幅スペクトルをプロット
plt.subplot(3, 1, 3)
plt.plot(f, P1_clean, linewidth=2)
plt.title("Single-Sided Amplitude Spectrum of Clean Signal")
plt.xlabel("Frequency [Hz]")
plt.ylabel("|P1(f)|")
plt.grid()

# プロットの表示
plt.tight_layout()
plt.show()


import numpy as np
import matplotlib.pyplot as plt

# 1. 時間領域の信号とフーリエ変換
# サンプリング周波数と時間ベクトルの設定
Ts = 1 / 50          # サンプリング周期
t = np.arange(0, 10, Ts)  # 0秒から10秒まで、サンプリング周期で分割された時間ベクトル

# 15 Hzと20 Hzの正弦波信号の生成
x = np.sin(2 * np.pi * 15 * t) + np.sin(2 * np.pi * 20 * t)

# 時間領域信号のプロット
plt.figure(figsize=(14, 12))
plt.subplot(3, 2, 1)
plt.plot(t, x)
plt.xlabel('Time (seconds)')
plt.ylabel('Amplitude')
plt.title('Time Domain Signal (15 Hz and 20 Hz)')
plt.grid()

# フーリエ変換の実行
y = np.fft.fft(x)
fs = 1 / Ts  # サンプリング周波数
n = len(y)  # 信号の長さ
f = np.arange(0, n) * fs / n  # 周波数ベクトル

# 振幅スペクトルのプロット
plt.subplot(3, 2, 2)
plt.plot(f, np.abs(y))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('Amplitude Spectrum of Signal')
plt.grid()

# FFTシフトを用いたスペクトルのプロット
yshift = np.fft.fftshift(y)
fshift = np.arange(-n/2, n/2) * fs / n  # 周波数ベクトル（シフト後）

plt.subplot(3, 2, 3)
plt.plot(fshift, np.abs(yshift))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('Shifted Amplitude Spectrum')
plt.grid()

# 位相スペクトルのプロット
theta = np.angle(yshift)
plt.subplot(3, 2, 4)
plt.plot(fshift, theta)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase (radians)')
plt.title('Phase Spectrum of Signal')
plt.grid()

# 2. ノイズを含む信号のパワースペクトルプロット
# ノイズを含む信号の生成
np.random.seed(0)  # 再現性のために乱数のシードを固定
xnoise = x + 2.5 * np.random.randn(len(t))

# ノイズを含む信号のフーリエ変換
ynoise = np.fft.fft(xnoise)
ynoiseshift = np.fft.fftshift(ynoise)

# パワースペクトルの計算
power = np.abs(ynoiseshift) ** 2 / n

# ノイズを含む信号のパワースペクトルのプロット
plt.subplot(3, 2, 5)
plt.plot(fshift, power)
plt.title('Power Spectrum of Noisy Signal')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power')
plt.grid()

# 3. 位相スペクトルの抽出
# 位相の計算とプロット（振幅が小さい成分を0にする）
tol = 1e-6
z = yshift.copy()
z[np.abs(z) < tol] = 0  # 振幅の小さい成分を0にする

theta = np.angle(z)

plt.subplot(3, 2, 6)
plt.plot(fshift, theta / np.pi)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase / pi")
plt.title("Phase Spectrum with Small Amplitude Set to Zero")
plt.grid()

plt.tight_layout()
plt.show()

# 4. 位相付き正弦波信号のフーリエ変換
# 15 Hzと40 Hzの位相付き信号の生成
fs_phase = 100
t_phase = np.arange(0, 1, 1/fs_phase)  # 1秒間の時間ベクトル
x_phase = np.cos(2 * np.pi * 15 * t_phase - np.pi / 4) + np.sin(2 * np.pi * 40 * t_phase)

# フーリエ変換
y_phase = np.fft.fft(x_phase)
z_phase = np.fft.fftshift(y_phase)

# 周波数ベクトルの定義
ly = len(y_phase)
f_phase = (-ly/2 + np.arange(0, ly)) / ly * fs_phase

# 振幅スペクトルと位相スペクトルのプロット
plt.figure(figsize=(12, 8))

# 振幅スペクトルのプロット
plt.subplot(2, 1, 1)
plt.stem(f_phase, np.abs(z_phase), use_line_collection=True)
plt.xlabel("Frequency (Hz)")
plt.ylabel("|y|")
plt.title("Double-Sided Amplitude Spectrum of x(t)")
plt.grid()

# 位相スペクトルのプロット
z_phase[np.abs(z_phase) < tol] = 0
theta_phase = np.angle(z_phase)

plt.subplot(2, 1, 2)
plt.stem(f_phase, theta_phase / np.pi, use_line_collection=True)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Phase / π")
plt.title("Phase Spectrum of x(t)")
plt.grid()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 1. パラメータ設定
Fspan = 8e3  # 分析周波数レンジ (例: 8 kHz)
N = 256  # サンプル点数 (例: 2のべき乗)
fs = Fspan * 2.56  # サンプリング周波数
T = N / fs  # 取込時間（時間長）
dt = 1 / fs  # サンプリングの時間間隔
L = N  # 分析ライン数

# 2. 矩形波の生成
t = np.linspace(0, T, N, endpoint=False)  # 時間ベクトル
frequency = 50  # 矩形波の周波数 (例: 50 Hz)
signal = np.sign(np.sin(2 * np.pi * frequency * t))  # 矩形波生成

# 3. FFT計算
Y = np.fft.fft(signal)  # フーリエ変換
Y_magnitude = np.abs(Y) / N  # 振幅スペクトルを計算
Y_magnitude = Y_magnitude[:N//2] * 2  # 正の周波数成分のみ使用

# 周波数軸の生成
freqs = np.fft.fftfreq(N, d=dt)[:N//2]

# 4. 時間領域信号のプロット
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(t * 1e3, signal, linewidth=2)  # 時間軸をms単位に変換してプロット
plt.title(f'Time Domain Signal (Rectangular Wave, Frequency = {frequency} Hz)')
plt.xlabel('Time [ms]')
plt.ylabel('Amplitude')
plt.grid()

# 5. 周波数スペクトルのプロット
plt.subplot(2, 1, 2)
plt.plot(freqs / 1e3, 20 * np.log10(Y_magnitude), linewidth=2)
plt.title('Power Spectrum')
plt.xlabel('Frequency [kHz]')
plt.ylabel('Magnitude [dB]')
plt.grid()

plt.tight_layout()
plt.show()

import cmath
import math

# 複素数を表すクラスの代わりに Python の complex 型を使用

def fft_time(x, n, idft):
    """
    FFT または IFFT を実行する関数。
    
    Args:
        x (list of complex): 入力データの複素数配列
        n (int): データの個数
        idft (int): 0 のとき FFT, 1 のとき IFFT
    
    Returns:
        list of complex: FFT または IFFT の結果
    """
    # データのビット反転順序を変更する
    j = 0
    for i in range(1, n):
        bit = n // 2
        while j >= bit:
            j -= bit
            bit //= 2
        j += bit
        if i < j:
            x[i], x[j] = x[j], x[i]

    # FFT のバタフライ演算
    m = 2
    while m <= n:
        step = -2j * cmath.pi / m
        if idft == 1:
            step = -step
        w_m = cmath.exp(step)
        for k in range(0, n, m):
            w = 1
            for j in range(m // 2):
                t = w * x[k + j + m // 2]
                u = x[k + j]
                x[k + j] = u + t
                x[k + j + m // 2] = u - t
                w *= w_m
        m *= 2

    # IFFT の場合は正規化
    if idft == 1:
        for i in range(n):
            x[i] /= n
    
    return x

def main():
    # サンプルデータ：複素数のリストを定義
    # 本来はファイルから読み込むデータを使用
    n = 8  # データの個数（2^M 形式）
    M = int(math.log2(n))
    
    # 入力データ（例として sin 波を使用）
    x = [complex(math.sin(2 * math.pi * i / n), 0) for i in range(n)]

    print("入力データ:")
    for i in range(n):
        print(f"x[{i}] = {x[i]}")
    
    # FFT を実行
    fft_result = fft_time(x[:], n, 0)
    
    print("\nFFT 結果:")
    for i in range(n):
        print(f"X[{i}] = {fft_result[i]}")

    # IFFT を実行
    ifft_result = fft_time(fft_result[:], n, 1)
    
    print("\nIFFT 結果（元のデータに戻る）:")
    for i in range(n):
        print(f"x[{i}] = {ifft_result[i]}")

if __name__ == "__main__":
    main()

import numpy as np
import matplotlib.pyplot as plt

# 周波数の定義
f1 = 697
f2 = 1209
fs = 8192  # サンプリング周波数の定義

# サンプル数の定義
n = np.arange(0, 0.1, 1/fs)

# 合成信号の生成
x = np.sin(2 * np.pi * f1 * n) + np.sin(2 * np.pi * f2 * n)

# 信号のフーリエ変換の計算
X = np.fft.fft(x)

# 振幅特性の計算
magnitude = np.abs(X)

# 信号の逆フーリエ変換の計算
xx = np.real(np.fft.ifft(X))

# 元信号のプロット用時間軸の設定
t = np.linspace(0, 0.01, 82)

# プロットの設定
plt.figure(figsize=(10, 8))

# 元信号のプロット
plt.subplot(3, 1, 1)
plt.plot(t, x[:len(t)])
plt.grid()
plt.xlabel('Time [sec]', fontsize=8)
plt.ylabel('Amplitude', fontsize=8)
plt.title('Original Signal, FFT Signal and IFFT Signal', fontsize=10)

# フーリエ変換した信号のプロット
frequencies = np.fft.fftfreq(len(X), d=1/fs)
plt.subplot(3, 1, 2)
plt.plot(frequencies, magnitude)
plt.axis([0, 1500, 0, 450])
plt.grid()
plt.xlabel('Frequency [Hz]', fontsize=8)
plt.ylabel('Magnitude', fontsize=8)

# 逆フーリエ変換した信号のプロット
plt.subplot(3, 1, 3)
plt.plot(t, xx[:len(t)], 'r-')
plt.grid()
plt.xlabel('Time [sec]', fontsize=8)
plt.ylabel('Amplitude', fontsize=8)

# プロットの表示
plt.tight_layout()
plt.show()

import numpy as np

# 複素数の型を扱うためのクラス定義
class ComplexDbl:
    def __init__(self, real=0.0, imag=0.0):
        self.real = real
        self.imag = imag

    def __add__(self, other):
        return ComplexDbl(self.real + other.real, self.imag + other.imag)

    def __sub__(self, other):
        return ComplexDbl(self.real - other.real, self.imag - other.imag)

    def __mul__(self, other):
        real = self.real * other.real - self.imag * other.imag
        imag = self.real * other.imag + self.imag * other.real
        return ComplexDbl(real, imag)

    def __str__(self):
        return f"({self.real:.4f}, {self.imag:.4f}j)"

# FFTを行う関数
def FFTSub(NW, x, cst):
    # ビット反転による並べ替え
    n = len(x)
    j = 0
    for i in range(1, n):
        bit = n >> 1
        while j >= bit:
            j -= bit
            bit >>= 1
        j += bit
        if i < j:
            x[i], x[j] = x[j], x[i]
    
    # FFTのバタフライ演算
    stage = 1
    while stage < n:
        step = stage << 1
        for k in range(0, stage):
            angle = -2j * np.pi * k / step
            W = np.exp(angle)  # 回転因子
            for i in range(k, n, step):
                j = i + stage
                temp = W * complex(x[j].real, x[j].imag)
                x[j] = ComplexDbl(x[i].real - temp.real, x[i].imag - temp.imag)
                x[i] = ComplexDbl(x[i].real + temp.real, x[i].imag + temp.imag)
        stage = step

# 逆FFTを行う関数
def FFTInvSub(NW, x, cst):
    # xの虚部の符号を反転
    for i in range(len(x)):
        x[i].imag = -x[i].imag
    
    # FFTを再度適用
    FFTSub(NW, x, cst)
    
    # 結果を1/NWでスケーリング
    r = 1.0 / NW
    for i in range(len(x)):
        x[i].real *= r
        x[i].imag *= r
    
    # 虚部の符号を元に戻す
    for i in range(len(x)):
        x[i].imag = -x[i].imag

# サンプル実行
NW = 8  # FFTのデータ点数（2の累乗）
x = [ComplexDbl(i, 0) for i in range(NW)]  # 入力データ（実数部分に値を設定）
cst = np.cos(2 * np.pi * np.arange(NW // 4) / NW)  # コサインテーブルの作成

# FFTの実行
print("=== Original Data ===")
for val in x:
    print(val)

FFTSub(NW, x, cst)

print("=== FFT Result ===")
for val in x:
    print(val)

# 逆FFTの実行
FFTInvSub(NW, x, cst)

print("=== Inverse FFT Result ===")
for val in x:
    print(val)

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, ifft
import librosa

# 1. ダミー音声信号の生成（母音を模したサイン波）
# サンプリングレートと時間軸の設定
fs = 16000  # サンプリングレート（16kHz）
duration = 1.0  # 秒
t = np.linspace(0, duration, int(fs * duration), endpoint=False)

# 周波数成分（例：500Hz、1000Hz、1500Hzのサイン波を混合して母音を模倣）
freqs = [500, 1000, 1500]
y = sum([np.sin(2 * np.pi * f * t) for f in freqs])

# 2. ハニング窓の適用
window = np.hanning(len(y))
y_windowed = y * window

# 3. フーリエ変換
fft_result = fft(y_windowed)
power_spectrum = np.abs(fft_result)**2  # パワースペクトル
log_power_spectrum = np.log(power_spectrum + 1e-10)  # 対数パワースペクトル

# 4. ケプストラムの計算
cepstrum = ifft(log_power_spectrum).real

# 5. 高次元成分の除去（リフタリング）
lifter = np.ones(len(cepstrum))
lifter[50:] = 0  # 低次元成分のみを残す
cepstrum_liftered = cepstrum * lifter

# 6. スペクトル包絡の抽出
envelope = fft(cepstrum_liftered).real

# 7. 平均スペクトル包絡のプロット
plt.figure(figsize=(10, 6))

# オリジナルのスペクトル
plt.subplot(2, 1, 1)
plt.title("Original Power Spectrum")
plt.plot(np.abs(fft_result))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")

# スペクトル包絡
plt.subplot(2, 1, 2)
plt.title("Spectral Envelope (Cepstrum Method)")
plt.plot(envelope)
plt.xlabel("Quefrency")
plt.ylabel("Amplitude")

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 1. ダミー信号の生成（矩形波）
fs = 1000  # サンプリングレート（1kHz）
T = 1.0  # 信号の周期（1秒）
t = np.linspace(0, T, int(fs * T), endpoint=False)  # 時間軸
signal = np.where((t % T) < 0.3, 1.0, 0.0)  # パルス幅0.3秒の矩形波

# 2. FFT の実行
N = 1024  # FFT のサンプル数（2のべき乗）
signal_padded = np.pad(signal, (0, N - len(signal)), 'constant')  # ゼロパディング

# フーリエ変換
fft_result = np.fft.fft(signal_padded)
frequencies = np.fft.fftfreq(N, 1/fs)

# 3. iFFT の実行（逆フーリエ変換）
ifft_result = np.fft.ifft(fft_result).real

# 4. 結果のプロット
plt.figure(figsize=(12, 8))

# オリジナルの信号
plt.subplot(3, 1, 1)
plt.plot(t, signal, label='Original Signal')
plt.title('Original Signal (Time Domain)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()

# フーリエ変換の結果
plt.subplot(3, 1, 2)
plt.plot(frequencies[:N//2], np.abs(fft_result)[:N//2], label='FFT (Magnitude)')
plt.title('FFT Result (Frequency Domain)')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude')
plt.legend()

# 逆フーリエ変換の結果
plt.subplot(3, 1, 3)
plt.plot(np.linspace(0, T, len(ifft_result)), ifft_result, label='Reconstructed Signal')
plt.title('Reconstructed Signal (iFFT)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.signal import get_window

# 1. ダミー信号の生成（2つの正弦波の加算）
fs = 1000  # サンプリングレート (Hz)
T = 1.0    # 信号の持続時間 (秒)
t = np.linspace(0, T, int(fs * T), endpoint=False)  # 時間軸
# 50Hzと120Hzの正弦波を合成
signal = 1.0 * np.sin(2.0 * np.pi * 50.0 * t) + 0.5 * np.sin(2.0 * np.pi * 120.0 * t)

# 2. 窓処理（ハニング窓、ハミング窓、ブラックマン窓）
windows = ['hann', 'hamming', 'blackman']
fig, axs = plt.subplots(len(windows), 2, figsize=(12, 8))

for i, window_name in enumerate(windows):
    # 窓関数を取得して信号に適用
    window = get_window(window_name, len(signal))
    windowed_signal = signal * window

    # FFTの計算
    fft_result = fft(windowed_signal)
    freqs = np.fft.fftfreq(len(fft_result), 1/fs)

    # 結果のプロット（時間領域）
    axs[i, 0].plot(t, windowed_signal)
    axs[i, 0].set_title(f'{window_name.capitalize()} Window - Time Domain')
    axs[i, 0].set_xlabel('Time [s]')
    axs[i, 0].set_ylabel('Amplitude')

    # 結果のプロット（周波数領域）
    axs[i, 1].plot(freqs[:len(freqs)//2], np.abs(fft_result)[:len(fft_result)//2])
    axs[i, 1].set_title(f'{window_name.capitalize()} Window - Frequency Domain')
    axs[i, 1].set_xlabel('Frequency [Hz]')
    axs[i, 1].set_ylabel('Magnitude')

plt.tight_layout()
plt.show()



import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.signal import get_window

# 1. アナログ信号（正弦波）を生成
fs = 1000  # 基本サンプリング周波数（Hz）
t = np.linspace(0, 1, fs, endpoint=False)  # 1秒間の時間軸
f = 5  # 信号の周波数（Hz）
analog_signal = np.sin(2 * np.pi * f * t)  # 5Hzの正弦波を生成

# 2. サンプリングと量子化のシミュレーション
quantization_levels = 16  # 量子化レベル（例: 16レベル = 4ビット）
quantized_signal = np.round((analog_signal + 1) * (quantization_levels / 2 - 1)) / (quantization_levels / 2 - 1) - 1

# 3. 窓関数の適用とFFTの計算
window_name = 'hann'
window = get_window(window_name, len(quantized_signal))
windowed_signal = quantized_signal * window

fft_result_analog = fft(analog_signal)
fft_result_quantized = fft(quantized_signal)
fft_result_windowed = fft(windowed_signal)
freqs = np.fft.fftfreq(len(fft_result_analog), 1/fs)

# 4. オーバーサンプリングのシミュレーション
oversampling_factor = 4  # オーバーサンプリング係数
oversampled_fs = fs * oversampling_factor
oversampled_t = np.linspace(0, 1, oversampled_fs, endpoint=False)
oversampled_signal = np.sin(2 * np.pi * f * oversampled_t)
oversampled_fft_result = fft(oversampled_signal)
oversampled_freqs = np.fft.fftfreq(len(oversampled_fft_result), 1/oversampled_fs)

# 5. サブサンプリング（ダウンサンプリング）のシミュレーション
downsampling_factor = 2  # サブサンプリング係数
downsampled_signal = oversampled_signal[::downsampling_factor]
downsampled_fs = oversampled_fs / downsampling_factor
downsampled_fft_result = fft(downsampled_signal)
downsampled_freqs = np.fft.fftfreq(len(downsampled_fft_result), 1/downsampled_fs)

# 6. 結果のプロット（時間領域と周波数領域）
fig, axs = plt.subplots(4, 2, figsize=(14, 14))

# オリジナルのアナログ信号（時間領域）
axs[0, 0].plot(t, analog_signal, label='Analog Signal')
axs[0, 0].set_title('Original Analog Signal (Time Domain)')
axs[0, 0].set_xlabel('Time [s]')
axs[0, 0].set_ylabel('Amplitude')
axs[0, 0].legend()

# オリジナル信号のFFT（周波数領域）
axs[0, 1].plot(freqs[:fs//2], np.abs(fft_result_analog)[:fs//2])
axs[0, 1].set_title('FFT of Original Analog Signal (Frequency Domain)')
axs[0, 1].set_xlabel('Frequency [Hz]')
axs[0, 1].set_ylabel('Magnitude')

# 量子化信号（時間領域）
axs[1, 0].plot(t, quantized_signal, label='Quantized Signal')
axs[1, 0].set_title(f'Quantized Signal ({quantization_levels}-levels)')
axs[1, 0].set_xlabel('Time [s]')
axs[1, 0].set_ylabel('Amplitude')
axs[1, 0].legend()

# 量子化信号のFFT（周波数領域）
axs[1, 1].plot(freqs[:fs//2], np.abs(fft_result_quantized)[:fs//2])
axs[1, 1].set_title(f'FFT of Quantized Signal ({quantization_levels}-levels)')
axs[1, 1].set_xlabel('Frequency [Hz]')
axs[1, 1].set_ylabel('Magnitude')

# 窓処理後の信号（時間領域）
axs[2, 0].plot(t, windowed_signal, label=f'Windowed Signal ({window_name} window)')
axs[2, 0].set_title(f'Windowed Signal ({window_name} window)')
axs[2, 0].set_xlabel('Time [s]')
axs[2, 0].set_ylabel('Amplitude')
axs[2, 0].legend()

# 窓処理後の信号のFFT（周波数領域）
axs[2, 1].plot(freqs[:fs//2], np.abs(fft_result_windowed)[:fs//2])
axs[2, 1].set_title(f'FFT of Windowed Signal ({window_name} window)')
axs[2, 1].set_xlabel('Frequency [Hz]')
axs[2, 1].set_ylabel('Magnitude')

# オーバーサンプリング後の信号のFFT（周波数領域）
axs[3, 0].plot(oversampled_freqs[:oversampled_fs//2], np.abs(oversampled_fft_result)[:oversampled_fs//2])
axs[3, 0].set_title(f'FFT of Oversampled Signal ({oversampling_factor}x)')
axs[3, 0].set_xlabel('Frequency [Hz]')
axs[3, 0].set_ylabel('Magnitude')

# サブサンプリング後の信号のFFT（周波数領域）
axs[3, 1].plot(downsampled_freqs[:int(downsampled_fs)//2], np.abs(downsampled_fft_result)[:int(downsampled_fs)//2])
axs[3, 1].set_title(f'FFT of Downsampled Signal (Downsampled by {downsampling_factor})')
axs[3, 1].set_xlabel('Frequency [Hz]')
axs[3, 1].set_ylabel('Magnitude')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.signal import get_window

# 1. Generating the Analog Signal (Sine Wave)
sampling_rate = 1000  # Sampling rate in Hz
time = np.linspace(0, 1, sampling_rate, endpoint=False)  # 1-second time frame
frequency = 5  # Frequency of the analog signal (Hz)
analog_signal = np.sin(2 * np.pi * frequency * time)  # Generating a 5Hz sine wave

# 2. Simulation of Sampling and Quantization
quantization_levels = 16  # Quantization levels (e.g., 16 levels = 4-bit resolution)
quantized_signal = np.round((analog_signal + 1) * (quantization_levels / 2 - 1)) / (quantization_levels / 2 - 1) - 1

# 3. Apply Window Function and Perform FFT
window_name = 'hann'
window = get_window(window_name, len(quantized_signal))
windowed_signal = quantized_signal * window

# Calculate FFT
fft_result_analog = fft(analog_signal)
fft_result_quantized = fft(quantized_signal)
fft_result_windowed = fft(windowed_signal)
frequencies = np.fft.fftfreq(len(fft_result_analog), 1/sampling_rate)

# 4. Plotting the Results (English Labels and Titles)
fig, axs = plt.subplots(3, 2, figsize=(14, 12))

# Original Analog Signal in Time Domain
axs[0, 0].plot(time, analog_signal, label='Original Analog Signal')
axs[0, 0].set_title('Original Analog Signal (Time Domain)')
axs[0, 0].set_xlabel('Time [s]')
axs[0, 0].set_ylabel('Amplitude')
axs[0, 0].legend()

# FFT of Original Analog Signal
axs[0, 1].plot(frequencies[:sampling_rate//2], np.abs(fft_result_analog)[:sampling_rate//2])
axs[0, 1].set_title('FFT of Original Analog Signal')
axs[0, 1].set_xlabel('Frequency [Hz]')
axs[0, 1].set_ylabel('Magnitude')

# Quantized Signal in Time Domain
axs[1, 0].plot(time, quantized_signal, label=f'Quantized Signal ({quantization_levels} levels)')
axs[1, 0].set_title('Quantized Signal (Time Domain)')
axs[1, 0].set_xlabel('Time [s]')
axs[1, 0].set_ylabel('Amplitude')
axs[1, 0].legend()

# FFT of Quantized Signal
axs[1, 1].plot(frequencies[:sampling_rate//2], np.abs(fft_result_quantized)[:sampling_rate//2])
axs[1, 1].set_title(f'FFT of Quantized Signal ({quantization_levels} levels)')
axs[1, 1].set_xlabel('Frequency [Hz]')
axs[1, 1].set_ylabel('Magnitude')

# Windowed Signal in Time Domain
axs[2, 0].plot(time, windowed_signal, label=f'Windowed Signal ({window_name} Window)')
axs[2, 0].set_title(f'Windowed Signal ({window_name} Window)')
axs[2, 0].set_xlabel('Time [s]')
axs[2, 0].set_ylabel('Amplitude')
axs[2, 0].legend()

# FFT of Windowed Signal
axs[2, 1].plot(frequencies[:sampling_rate//2], np.abs(fft_result_windowed)[:sampling_rate//2])
axs[2, 1].set_title(f'FFT of Windowed Signal ({window_name} Window)')
axs[2, 1].set_xlabel('Frequency [Hz]')
axs[2, 1].set_ylabel('Magnitude')

# Display all plots
plt.tight_layout()
plt.show()