「三角関数和の周波数解析：基音440Hzに基づくFFT検証

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks

# ================= Parameters =================
fs = 44100        # Sampling frequency [Hz]
T = 1.0           # Signal duration [s] -> longer for higher frequency resolution
n = 10            # Number of harmonics
f0 = 440          # Base frequency (Hz, A4)

# ================= Time Axis =================
t = np.linspace(0, T, int(fs*T), endpoint=False)
sin_sum = np.zeros_like(t)
cos_sum = np.zeros_like(t)

# ================= Summation =================
for k in range(1, n+1):
    sin_sum += np.sin(2*np.pi*k*f0*t)
    cos_sum += np.cos(2*np.pi*k*f0*t)

# ================= Time-Domain Plot =================
plt.figure(figsize=(10,5))
plt.plot(t[:1000], sin_sum[:1000], label=r'$\sum_{k=1}^n \sin(2\pi k f_0 t)$')  # zoom in
plt.plot(t[:1000], cos_sum[:1000], label=r'$\sum_{k=1}^n \cos(2\pi k f_0 t)$')
plt.title(f"Trigonometric Harmonic Sums (n={n}, f0={f0} Hz)")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.show()

# ================= FFT Analysis (Zero Padding) =================
N = 2**15  # FFT size with zero padding for high resolution
sin_fft = np.fft.fft(sin_sum, n=N)
cos_fft = np.fft.fft(cos_sum, n=N)
freqs = np.fft.fftfreq(N, 1/fs)

# Keep only positive frequencies
pos_mask = freqs >= 0
freqs = freqs[pos_mask]
sin_fft_mag = np.abs(sin_fft)[pos_mask]
cos_fft_mag = np.abs(cos_fft)[pos_mask]

# ================= Peak Detection =================
threshold = 0.1  # relative threshold for peak detection
sin_peaks, _ = find_peaks(sin_fft_mag, height=np.max(sin_fft_mag)*threshold)
cos_peaks, _ = find_peaks(cos_fft_mag, height=np.max(cos_fft_mag)*threshold)
detected_sin_hz = freqs[sin_peaks]
detected_cos_hz = freqs[cos_peaks]

# ================= Frequency to Note Conversion =================
def freq_to_note(freq):
    if freq <= 0:
        return "N/A"
    midi_num = 69 + 12*np.log2(freq/440.0)
    midi_int = int(round(midi_num))
    note_names = ['C', 'C#', 'D', 'D#', 'E', 'F', 'F#', 'G', 'G#', 'A', 'A#', 'B']
    note = note_names[midi_int % 12] + str(midi_int//12 - 1)
    return note

# ================= FFT Plot with Peaks =================
plt.figure(figsize=(10,5))
plt.plot(freqs, sin_fft_mag, label="FFT of Σ sin(k·2πf0t)")
plt.plot(freqs, cos_fft_mag, label="FFT of Σ cos(k·2πf0t)")
plt.scatter(detected_sin_hz, sin_fft_mag[sin_peaks], color='red', marker='o', label="Detected peaks (sin)")
plt.scatter(detected_cos_hz, cos_fft_mag[cos_peaks], color='green', marker='x', label="Detected peaks (cos)")
plt.title(f"FFT Spectrum with Detected Harmonics (n={n}, f0={f0} Hz)")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Magnitude")
plt.xlim(0, f0*(n+1))  # show only up to highest harmonic
plt.legend()
plt.grid(True)
plt.show()

# ================= Print Detected Frequencies =================
print("Detected frequencies in Σsin(k·2πf0t):")
for f in detected_sin_hz:
    print(f"  {f:.2f} Hz -> {freq_to_note(f)}")

print("\nDetected frequencies in Σcos(k·2πf0t):")
for f in detected_cos_hz:
    print(f"  {f:.2f} Hz -> {freq_to_note(f)}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from IPython.display import Audio

# ================= Parameters =================
fs = 44100        # サンプリング周波数 [Hz]
T = 2.0           # 信号継続時間 [秒]
n = 10            # 倍音の数
f0 = 440          # 基音 A4 = 440Hz

# ================= 波形生成 =================
t = np.linspace(0, T, int(fs*T), endpoint=False)
sin_sum = np.zeros_like(t)
for k in range(1, n+1):
    sin_sum += np.sin(2*np.pi*k*f0*t)

# 正規化（音割れ防止）
sin_sum = sin_sum / np.max(np.abs(sin_sum))

# ================= 再生用埋め込みオーディオ =================
audio = Audio(sin_sum, rate=fs)
display(audio)

# ================= プロット =================
plt.figure(figsize=(10,4))
plt.plot(t[:2000], sin_sum[:2000])
plt.title(f"Harmonic Sum Signal (n={n}, f0={f0} Hz)")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()