FFTコード AD変換 Pythonコードメモ

import numpy as np
import matplotlib.pyplot as plt

# Constants and parameters
FFT_points = 512
fs = 3 * 10**6
Vth = 0.5
adc_bit = 8  # 8-bit ADC simulation
Radix = 2

Ts = 1 / fs
nyquist_points = FFT_points // 2
t = np.arange(0, FFT_points * Ts, Ts)
s = np.linspace(1, nyquist_points, nyquist_points - 1)

# Simulate 8-bit ADC output
ad_output = np.random.randint(0, 256, size=FFT_points)  # Simulated ADC output (8-bit)
ad_signal = (ad_output / 255.0) * (2 * Vth) - Vth  # Normalize ADC output to original signal range

# Simulate DAC conversion with added noise
noise_amplitude = 0.1  # Amplitude of simulated noise
da_signal = ad_signal + np.random.normal(0, noise_amplitude, FFT_points)  # Add noise to simulate DAC output

# FFT processing
xfft_da = np.fft.rfft(da_signal) * 2 / FFT_points
P_da = (xfft_da * xfft_da.conj()).real
P_dB_da = 10 * np.log10(P_da)

# Calculate SNDR and ENOB based on DAC output
fin_index_da = np.argmax(P_da[:nyquist_points])
SNDR_da = 10 * np.log10(P_da[fin_index_da] / (sum(P_da[:nyquist_points]) - P_da[fin_index_da]))
ENOB_da = (SNDR_da - 1.76) / 6.02

# Plotting
fig4 = plt.figure()
plt.title("DAC Conversion Output (" + str(FFT_points) + " points)")
plt.xlabel("Time [s]")
plt.ylabel("Voltage [V]")
plt.grid(True)
plt.plot(t, da_signal, drawstyle='steps-post')

fig5 = plt.figure()
plt.title("Power Spectrum of DAC Output (" + str(FFT_points) + " points)")
plt.xlabel("Normalized Frequency (f/fs)")
plt.ylabel("Power [dB]")
plt.text(0.3, -10, "SNDR = " + format(SNDR_da, '.2f') + "[dB]", size=10)
plt.text(0.3, -20, "ENOB = " + format(ENOB_da, '.2f') + "[bit]", size=10)
plt.grid(True)
plt.plot(s / FFT_points, P_dB_da[1:nyquist_points])

fig6 = plt.figure()
plt.title("Normalized DAC Output (" + str(FFT_points) + " points)")
plt.xlabel("Time [s]")
plt.ylabel("Code")
plt.grid(True)
plt.scatter(t, da_signal, s=1)

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft

# Parameters
Fs = 1000  # Sampling frequency (Hz)
T = 1 / Fs  # Sampling period
t = np.arange(0, 1, T)  # Time vector
f1 = 5  # Frequency of the sinusoid (Hz)
A = 1  # Amplitude of the sinusoid

# Generate sinusoidal signal
sin_wave = A * np.sin(2 * np.pi * f1 * t)

# Simulate ADC (quantization to 8 bits)
quantized_signal = np.round((sin_wave + 1) * 127.5) / 127.5
quantized_signal = np.clip(quantized_signal, -1, 1)  # Clip to valid range [-1, 1]

# Compute FFT
fft_signal = fft(quantized_signal)
power_spectrum = np.abs(fft_signal)**2

# Calculate SNR
signal_power = np.sum(np.abs(fft(sin_wave))**2)  # Power of the original signal
noise_power = np.sum(power_spectrum) - signal_power  # Power of the noise

SNR_dB = 10 * np.log10(signal_power / noise_power)

# Calculate frequency bins for plotting
N = len(t)
freq = np.fft.fftfreq(N, T)[:N//2]  # FFT frequency bins
power_spectrum_plot = power_spectrum[:N//2] / N  # Normalize by N and take only positive frequencies

# Plot
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(t, sin_wave, label='Original Signal')
plt.plot(t, quantized_signal, label='Quantized Signal')
plt.title('Original and Quantized Signals')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(freq, 10 * np.log10(power_spectrum_plot))
plt.title('Power Spectrum (FFT)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (dB)')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.text(0.1, 0.5, f'SNR: {SNR_dB:.2f} dB', fontsize=12, ha='center')
plt.axis('off')

plt.tight_layout()
plt.show()

print(f"SNR: {SNR_dB:.2f} dB")

import numpy as np
import matplotlib.pyplot as plt

# Parameters
fs = 80e6  # Sampling frequency (80.000 MSPS)
fiN = 2.111e6  # Input frequency (2.111 MHz)
duration = 1.0  # Duration for analysis (seconds)
t = np.linspace(0, duration, int(fs * duration), endpoint=False)  # Time axis

# Generate input signal (sine wave)
input_signal = np.sin(2 * np.pi * fiN * t)

# Quantization (12-bit ADC)
quantized_signal = np.round(input_signal * (2**11)) / (2**11)

# Compute FFT
fft_output = np.fft.fft(quantized_signal, 4096)
freq_axis = np.fft.fftfreq(4096, d=1/fs) / 1e6  # Frequency axis (MHz)

# Calculate Power Spectral Density (PSD)
psd = 20 * np.log10(np.abs(fft_output))

# Calculate SFDR and SNR
fundamental_idx = np.argmax(psd)  # Index of the fundamental frequency
psd[fundamental_idx] = -np.inf  # Ignore power of the fundamental frequency to compute SFDR
SFDR = np.max(psd)  # SFDR is the height of the highest sidelobe
SNR = psd[fundamental_idx]  # SNR is the power of the fundamental frequency

# Plot
plt.figure(figsize=(10, 6))

# Plot PSD
plt.plot(freq_axis, psd, label='PSD (dB)')
plt.axhline(y=SFDR, color='r', linestyle='--', label=f'SFDR = {SFDR:.1f} dBc')
plt.axhline(y=SNR, color='g', linestyle='--', label=f'SNR = {SNR:.1f} dBc')

plt.title('Analysis of Ideal 12-bit ADC Output with 4096-point FFT')
plt.xlabel('Frequency (MHz)')
plt.ylabel('Power Spectral Density (dB)')
plt.xlim(0, 40)
plt.grid(True)
plt.legend()
plt.tight_layout()

plt.show()

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

# サンプリング周波数と時間軸の設定
fs = 1000  # サンプリング周波数
t = np.arange(0, 1, 1/fs)  # 1秒間の時間軸
f_signal = 10  # サイン波の周波数
signal = np.sin(2 * np.pi * f_signal * t)  # サイン波信号

# 6ビットのDAC出力にノイズを加える（ここでは単純にサイン波にガウシアンノイズを加える例）
noise_power = 0.01 * np.var(signal)  # ノイズのパワー（信号の分散の1%）
noise = np.random.normal(0, np.sqrt(noise_power), signal.shape)  # ガウシアンノイズ
signal_noisy = signal + noise  # ノイズが乗ったサイン波信号

# FFTを計算
fft_signal = fft(signal_noisy)
freqs = fftfreq(len(signal_noisy), 1/fs)

# パワースペクトルを計算
power_spectrum = np.abs(fft_signal)**2

# メインローブのピークを見つける
peaks, _ = find_peaks(power_spectrum[:len(freqs)//2], height=100)  # 高さの閾値は適宜調整

# SNRを計算（ピークを除いた平均パワーとノイズパワーの比）
signal_power = np.sum(power_spectrum[peaks[0]])  # メインローブのピークのパワー
noise_power_fft = np.sum(power_spectrum) - signal_power  # ピーク以外の全ての周波数成分のパワー
SNR = 10 * np.log10(signal_power / noise_power_fft)

# SNDRを計算（サイン波パワーとノイズフロアの比）
noise_floor = np.mean(power_spectrum[len(freqs)//2:])
SNDR = 10 * np.log10(signal_power / noise_floor)

# SFDRを計算（最大ピークとノイズフロアの比）
max_peak = np.max(power_spectrum[peaks])
SFDR = 10 * np.log10(max_peak / noise_floor)

# ENOBを計算（SNDRから）
ENOB = (SNDR - 1.76) / 6.02

# 結果の出力
print(f"SNR: {SNR} dB")
print(f"SNDR: {SNDR} dB")
print(f"SFDR: {SFDR} dB")
print(f"ENOB: {ENOB}")

# FFTのプロット
plt.figure(figsize=(12, 6))
plt.plot(freqs[:len(freqs)//2], 10*np.log10(power_spectrum[:len(freqs)//2]), label='Power Spectrum')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power (dB)')
plt.title('FFT of Noisy 6-bit DAC Output Signal')
plt.legend()
plt.grid(True)
plt.show()

import matplotlib.pyplot as plt
import numpy as np

# Constants
FFT_points = 4096  # FFT points
fin = 10           # [Hz] Input frequency
Ampl = 1           # Input signal amplitude

# Sampling frequency calculation
fs = fin * (FFT_points) / (2**6 - 1)
print(f"Sampling Frequency: {fs} Hz")

# Define the sine wave function
def sin_wave(t):
    return Ampl * np.sin(2 * np.pi * fin * t)

# Generate the time axis for plotting
tp = np.linspace(0, 0.1, 1000)  # Time axis for plotting

# Plot the input signal
plt.figure()
plt.title("Input Signal")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.xlim(0, 0.1)  # Plot only the 0 to 0.1 second interval
plt.plot(tp, sin_wave(tp))
plt.show()

# ADC parameters
bit = 9
D = [1]  # Initialize MSB to 1
zeros = [0] * (bit - 1)  # Set other bits to 0
D.extend(zeros)  # Initialize register with the specified values

# DAC conversion
weight = 2**np.arange(0, -bit, -1, dtype=float)
dac_out = np.dot(D, weight)

# Generate the signal for FFT
t = np.arange(0, FFT_points) / fs  # Time axis for the signal
y = sin_wave(t)

# FFT calculation and normalization
F = np.fft.rfft(y) / FFT_points * 2  # FFT calculation and normalization

# Power calculation
P = (F * F.conj()).real  # Power calculation
with np.errstate(divide='ignore', invalid='ignore'):
    P_dB = 10 * np.log10(P)  # Convert power to dB
P_dB[P_dB == -np.inf] = -200  # Replace -inf with -200

# SNR calculation
fin_index = int(np.argmax(P[:FFT_points // 2]))
signal_power = P[fin_index]
noise_power = sum(P[:(FFT_points // 2)]) - signal_power
noise_power = max(noise_power, 1e-20)  # Add small constant to avoid division by zero
SNR = 10 * np.log10(signal_power / noise_power)
ENOB = (SNR - 1.76) / 6.02

# Print results
print(f"Bit array: {D}")
print(f"DAC Output: {dac_out}")
print(f"SNR: {SNR} dB")
print(f"ENOB: {ENOB}")

# Plot the FFT result
freqs = np.fft.rfftfreq(FFT_points, 1/fs)
plt.figure()
plt.plot(freqs, P_dB)
plt.title("FFT Result")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power [dB]")
plt.grid(True)
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
k = 1.38e-23  # ボルツマン定数 (J/K)
q = 1.6e-19   # 電子の電荷 (C)
T = 300       # 絶対温度 (K)
R = 1000      # 抵抗値 (Ω)
I = 1e-6      # 電流 (A)
alpha = 1.2   # 1/fノイズの指数

# 周波数の範囲
frequencies = np.logspace(1, 6, 100)

# 熱ノイズのパワー密度
thermal_noise_psd = 4 * k * T * R

# ショットノイズのパワー密度
shot_noise_psd = 2 * q * I

# 1/fノイズのパワー密度
flicker_noise_psd = 1 / frequencies**alpha

# プロット
plt.figure(figsize=(12, 8))

# 熱ノイズのプロット
plt.axhline(y=thermal_noise_psd, color='r', linestyle='--', label='Johnson Noise (Thermal)')

# ショットノイズのプロット
plt.axhline(y=shot_noise_psd, color='b', linestyle='--', label='Shot Noise')

# 1/fノイズのプロット
plt.loglog(frequencies, flicker_noise_psd, label='1/f Noise', color='g')

plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (PSD)')
plt.title('Noise Power Spectral Density')
plt.legend()
plt.grid(True, which="both", ls="--")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# 定数
k = 1.38e-23  # ボルツマン定数 (J/K)
T = 300       # 絶対温度 (K)
R = 1000      # 抵抗値 (Ω)
q = 1.6e-19   # 電子の電荷 (C)
I = 1e-6      # 電流 (A)
alpha = 1.2   # 1/fノイズの指数
signal_power = 1e-12  # 信号のパワー (W)

# 周波数帯域幅
bandwidth = np.logspace(1, 6, 10)  # 10 Hz から 1 MHz までの帯域幅

# ノイズパワー密度の計算
thermal_noise_psd = 4 * k * T * R
shot_noise_psd = 2 * q * I
frequencies = np.logspace(1, 6, 100)
flicker_noise_psd = 1 / frequencies**alpha

# ノイズパワーの計算
thermal_noise_power = thermal_noise_psd * bandwidth
shot_noise_power = shot_noise_psd * bandwidth

# SNRの計算
snr_thermal = 10 * np.log10(signal_power / thermal_noise_power)
snr_shot = 10 * np.log10(signal_power / shot_noise_power)

# プロット
plt.figure(figsize=(14, 6))

# ノイズパワーのプロット
plt.subplot(1, 2, 1)
plt.loglog(bandwidth, thermal_noise_power, label='Thermal Noise Power')
plt.loglog(bandwidth, shot_noise_power, label='Shot Noise Power')
plt.xlabel('Bandwidth (Hz)')
plt.ylabel('Noise Power (W)')
plt.title('Noise Power vs Bandwidth')
plt.legend()
plt.grid(True, which="both", ls="--")

# SNRのプロット
plt.subplot(1, 2, 2)
plt.plot(bandwidth, snr_thermal, label='SNR (Thermal Noise)')
plt.plot(bandwidth, snr_shot, label='SNR (Shot Noise)')
plt.xlabel('Bandwidth (Hz)')
plt.ylabel('SNR (dB)')
plt.title('SNR vs Bandwidth')
plt.legend()
plt.grid(True, which="both", ls="--")

plt.tight_layout()
plt.show()