convert-circuitjs-export-text-to-wav.py with sinc1-5 filter

convert-circuitjs-export-text-to-wav.py

from scipy.io import wavfile
import scipy
import matplotlib.pyplot as plt
import numpy as np
import sys

def simple_moving_average_numpy(data, window_size):
    """Calculates the Simple Moving Average using numpy.convolve."""
    # Create an impulse response (weights) array where each weight is 1/window_size
    weights = np.ones(window_size) / window_size
    # Convolve the data with the weights using 'valid' mode to return only full overlaps
    sma = np.convolve(data, weights, mode='valid')
    return sma

print("reading circuitjs.")


filename_in = "/home/e/Downloads/data-20251211-1914.circuitjs.txt" #sys.argv[1]

filelines = open(filename_in).readlines()

del filelines[0]


#data = np.int16(float(filelines)*65536/2/5)
data = np.zeros(len(filelines))
#np.empty([1], dtype=int)

maxint16 = np.iinfo(np.int16).max

for i in range(len(filelines)):
    value = float(filelines[i]);
    #print(i, value)
    data[i] = (value - 2.5)

print(data)
# Calculate the mean of the entire array
average = np.mean(data)
print(f"Average of the entire array: {average}")


# Specify the path to your WAV file
filename_original_out = filename_in + ".wav"

# Write the data to a WAV file
fs_original = 100000000
simple_moving_average_data1 = simple_moving_average_numpy(data, fs_original // 40000)
simple_moving_average_data2 = simple_moving_average_numpy(simple_moving_average_data1, fs_original // 40000)
simple_moving_average_data3 = simple_moving_average_numpy(simple_moving_average_data2, fs_original // 40000)
simple_moving_average_data4 = simple_moving_average_numpy(simple_moving_average_data3, fs_original // 40000)
simple_moving_average_data5 = simple_moving_average_numpy(simple_moving_average_data4, fs_original // 40000)
decimation_factor = 50
fs_decimated = int(fs_original / decimation_factor)
decimated_data0 = scipy.signal.decimate(data, q=decimation_factor, ftype='fir', zero_phase=True)
decimated_data1 = scipy.signal.decimate(simple_moving_average_data1, q=decimation_factor, ftype='fir', zero_phase=True)
decimated_data2 = scipy.signal.decimate(simple_moving_average_data2, q=decimation_factor, ftype='fir', zero_phase=True)
decimated_data3 = scipy.signal.decimate(simple_moving_average_data3, q=decimation_factor, ftype='fir', zero_phase=True)
decimated_data4 = scipy.signal.decimate(simple_moving_average_data4, q=decimation_factor, ftype='fir', zero_phase=True)
decimated_data5 = scipy.signal.decimate(simple_moving_average_data5, q=decimation_factor, ftype='fir', zero_phase=True)

wavfile.write(filename_original_out, fs_decimated, decimated_data0)
print("WAV file created successfully.")

plt.figure()
plt.psd(decimated_data0, Fs=fs_decimated, NFFT=2048, noverlap=128, scale_by_freq=True, label='sinc0')
plt.psd(decimated_data1, Fs=fs_decimated, NFFT=2048, noverlap=128, scale_by_freq=True, label='sinc1')
plt.psd(decimated_data2, Fs=fs_decimated, NFFT=2048, noverlap=128, scale_by_freq=True, label='sinc2')
plt.psd(decimated_data3, Fs=fs_decimated, NFFT=2048, noverlap=128, scale_by_freq=True, label='sinc3')
plt.psd(decimated_data4, Fs=fs_decimated, NFFT=2048, noverlap=128, scale_by_freq=True, label='sinc4')
plt.psd(decimated_data5, Fs=fs_decimated, NFFT=2048, noverlap=128, scale_by_freq=True, label='sinc5')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (dB/Hz)')
plt.ylim(-120, 0)
plt.yticks([0, -10, -20, -30, -40, -50, -60, -70, -80, -90, -100, -110, -120, -130, -140, -150, -160])
plt.title('Power Spectral Density')
plt.grid(True)
plt.legend()
plt.show()

'''
# Calculate the PSD using Welch's method
frequencies, psd = scipy.signal.welch(decimated_data, fs_decimated, window='blackmanharris', nperseg=8192)

# Plot the PSD
plt.figure()
plt.semilogy(frequencies, psd)
plt.title('Power Spectral Density')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD (V^2/Hz)')
plt.grid(True)
plt.show()
'''
'''
# f contains the frequency components
# S is the PSD
(f, S) = scipy.signal.periodogram(decimated_data, fs_decimated, scaling='density')

plt.semilogy(f, S)
plt.ylim([1e-9, 1e0])
plt.xlim([0,fs_decimated/2])
plt.xlabel('frequency [Hz]')
plt.ylabel('PSD [V**2/Hz]')
plt.show()

'''




'''
decimated_data = signal.decimate(data, q=100, ftype='fir', zero_phase=True)

# Plot the Power Spectral Density (PSD)
plt.figure(figsize=(10, 6))
plt.psd(decimated_data, NFFT=256, Fs=100000, Fc=0, scale_by_freq=True)
plt.title('Power Spectral Density (PSD)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power/Frequency (dB/Hz)')
plt.grid(True)
plt.show()


frequencies2, psd2 = signal.welch(decimated_data, 1000000, nperseg=262144)
Psd2_db = 10 * np.log10(psd2)

plt.figure()
plt.plot(frequencies2, Psd2_db) # semilogy sets y-axis to log by default
plt.xscale('log') # Set x-axis to log scale
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (dB/Hz)')
plt.title('PSD with Logarithmic Frequency Axis')
plt.grid(True, which="both", ls="-")
plt.show()

# Generate a sample signal (e.g., a sine wave with some noise)
fs = 100000000  # Sampling frequency (Hz)
t = np.arange(0, 1, 1/fs)  # Time vector from 0 to 1 second
#noise = np.random.randn(len(t)) * 0.2  # Add some random noise


frequencies, psd = signal.welch(data, fs, nperseg=1024)

# Define the selected frequency band (e.g., 20 Hz to 20 kHz)
low_freq_band = 20
high_freq_band = 20000

# Find the indices corresponding to the selected frequency band
idx_band = np.logical_and(frequencies >= low_freq_band, frequencies <= high_freq_band)

# Extract the PSD and frequencies within the selected band
freqs_band = frequencies[idx_band]
psd_band = psd[idx_band]

# Plot the full PSD and highlight the selected band
plt.figure(figsize=(10, 6))
#plt.semilogy(frequencies, psd, label='Full PSD')
plt.semilogy(freqs_band, psd_band, 'r', linewidth=2, label=f'Selected Band ({low_freq_band}-{high_freq_band} Hz)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (V^2/Hz)')
plt.title('Power Spectral Density with Selected Frequency Band')
plt.legend()
plt.grid(True)
plt.show()

plt.figure(figsize=(10, 6))
plt.semilogx(frequencies, psd) # Use semilogx for log-scaled x-axis
plt.title('Power Spectral Density with Logarithmic Frequency Axis')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power/Frequency (dB/Hz)')
plt.grid(True, which="both", ls="-") # Add grid lines for clarity
plt.show()

# Plot the Power Spectral Density (PSD)
plt.figure(figsize=(10, 6))
plt.psd(data, NFFT=256, Fs=fs, Fc=0, scale_by_freq=True)
plt.title('Power Spectral Density (PSD)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power/Frequency (dB/Hz)')
plt.grid(True)
plt.show()
'''

lower null to 40kHz:

looking just at low frequencies: