Export FIR coefficients for analysis

Author: kunitoki

tests/yup_dsp/yup_FilterDesigner.cpp

@@ -516,3 +516,152 @@ TEST_F (FilterDesignerTests, FirFloatDoubleConsistency)
     for (int i = 0; i < numCoeffs; ++i)
         EXPECT_NEAR (doubleCoeffs[i], static_cast<double> (floatCoeffs[i]), toleranceF);
 }
+
+TEST_F (FilterDesignerTests, DISABLED_ExportFIRCoefficientsForAnalysis)
+{
+    const int numCoeffs = 129;
+    const float sampleRateF = 44100.0f;
+
+    // Design different FIR filters
+    auto lowpass = FilterDesigner<float>::designFIRLowpass (numCoeffs, 10000.0f, sampleRateF);
+    auto highpass = FilterDesigner<float>::designFIRHighpass (numCoeffs, 10000.0f, sampleRateF);
+    auto bandpass = FilterDesigner<float>::designFIRBandpass (numCoeffs, 8000.0f, 12000.0f, sampleRateF);
+    auto bandstop = FilterDesigner<float>::designFIRBandstop (numCoeffs, 8000.0f, 12000.0f, sampleRateF);
+
+    // Different windows for lowpass
+    auto lowpassHann = FilterDesigner<float>::designFIRLowpass (numCoeffs, 10000.0f, sampleRateF, WindowType::hann);
+    auto lowpassHamming = FilterDesigner<float>::designFIRLowpass (numCoeffs, 10000.0f, sampleRateF, WindowType::hamming);
+    auto lowpassBlackman = FilterDesigner<float>::designFIRLowpass (numCoeffs, 10000.0f, sampleRateF, WindowType::blackman);
+
+    // Helper lambda to write coefficients to file
+    auto writeCoeffs = [] (const std::vector<float>& coeffs, const std::string& filename)
+    {
+        std::ofstream file (filename);
+        if (file.is_open())
+        {
+            for (size_t i = 0; i < coeffs.size(); ++i)
+            {
+                file << coeffs[i];
+                if (i < coeffs.size() - 1)
+                    file << "\n";
+            }
+            file.close();
+        }
+    };
+
+    // Write all coefficient sets to files
+    writeCoeffs (lowpass, "fir_lowpass_10000hz.txt");
+    writeCoeffs (highpass, "fir_highpass_10000hz.txt");
+    writeCoeffs (bandpass, "fir_bandpass_8000_12000hz.txt");
+    writeCoeffs (bandstop, "fir_bandstop_8000_12000hz.txt");
+    writeCoeffs (lowpassHann, "fir_lowpass_hann_10000hz.txt");
+    writeCoeffs (lowpassHamming, "fir_lowpass_hamming_10000hz.txt");
+    writeCoeffs (lowpassBlackman, "fir_lowpass_blackman_10000hz.txt");
+
+    // Create a Python script to plot the frequency responses
+    std::ofstream pyScript ("plot_fir_responses.py");
+    if (pyScript.is_open())
+    {
+        pyScript << R"(#!/usr/bin/env python3
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy import signal
+
+def load_coeffs(filename):
+    with open(filename, 'r') as f:
+        return [float(line.strip()) for line in f.readlines()]
+
+def plot_frequency_response(coeffs, title, sample_rate=44100):
+    w, h = signal.freqz(coeffs, worN=8000, fs=sample_rate)
+    
+    plt.figure(figsize=(12, 8))
+    
+    # Magnitude response
+    plt.subplot(2, 1, 1)
+    plt.plot(w, 20 * np.log10(np.abs(h)))
+    plt.title(f'{title} - Magnitude Response')
+    plt.xlabel('Frequency (Hz)')
+    plt.ylabel('Magnitude (dB)')
+    plt.grid(True)
+    plt.xlim(0, sample_rate/2)
+    plt.ylim(-80, 5)
+    
+    # Phase response
+    plt.subplot(2, 1, 2)
+    plt.plot(w, np.unwrap(np.angle(h)) * 180 / np.pi)
+    plt.title(f'{title} - Phase Response')
+    plt.xlabel('Frequency (Hz)')
+    plt.ylabel('Phase (degrees)')
+    plt.grid(True)
+    plt.xlim(0, sample_rate/2)
+    
+    plt.tight_layout()
+    plt.savefig(f'{title.lower().replace(" ", "_").replace("-", "_")}_response.png', dpi=150, bbox_inches='tight')
+    plt.show()
+
+# Load and plot all FIR filter responses
+filters = [
+    ('fir_lowpass_10000hz.txt', 'FIR Lowpass 10000Hz'),
+    ('fir_highpass_10000hz.txt', 'FIR Highpass 10000Hz'), 
+    ('fir_bandpass_8000_12000hz.txt', 'FIR Bandpass 8000-12000Hz'),
+    ('fir_bandstop_8000_12000hz.txt', 'FIR Bandstop 8000-12000Hz'),
+    ('fir_lowpass_hann_10000hz.txt', 'FIR Lowpass Hann Window'),
+    ('fir_lowpass_hamming_10000hz.txt', 'FIR Lowpass Hamming Window'),
+    ('fir_lowpass_blackman_10000hz.txt', 'FIR Lowpass Blackman Window')
+]
+
+for filename, title in filters:
+    try:
+        coeffs = load_coeffs(filename)
+        plot_frequency_response(coeffs, title)
+    except FileNotFoundError:
+        print(f"File {filename} not found!")
+
+# Compare window types on same plot
+plt.figure(figsize=(12, 6))
+window_files = [
+    ('fir_lowpass_hann_1000hz.txt', 'Hann', 'blue'),
+    ('fir_lowpass_hamming_1000hz.txt', 'Hamming', 'red'),
+    ('fir_lowpass_blackman_1000hz.txt', 'Blackman', 'green')
+]
+
+for filename, label, color in window_files:
+    try:
+        coeffs = load_coeffs(filename)
+        w, h = signal.freqz(coeffs, worN=8000, fs=44100)
+        plt.plot(w, 20 * np.log10(np.abs(h)), label=label, color=color)
+    except FileNotFoundError:
+        print(f"File {filename} not found!")
+
+plt.title('FIR Lowpass 1000Hz - Window Comparison')
+plt.xlabel('Frequency (Hz)')
+plt.ylabel('Magnitude (dB)')
+plt.grid(True)
+plt.legend()
+plt.xlim(0, 22050)
+plt.ylim(-80, 5)
+plt.savefig('fir_window_comparison.png', dpi=150, bbox_inches='tight')
+plt.show()
+
+print("All plots generated successfully!")
+)";
+        pyScript.close();
+    }
+
+    // Just verify the files were created - the actual validation will be done visually with Python
+    EXPECT_EQ (lowpass.size(), numCoeffs);
+    EXPECT_EQ (highpass.size(), numCoeffs);
+    EXPECT_EQ (bandpass.size(), numCoeffs);
+    EXPECT_EQ (bandstop.size(), numCoeffs);
+
+    std::cout << "\nFIR coefficient files and Python plotting script created:\n";
+    std::cout << "- fir_lowpass_10000hz.txt\n";
+    std::cout << "- fir_highpass_10000hz.txt\n";
+    std::cout << "- fir_bandpass_8000_12000hz.txt\n";
+    std::cout << "- fir_bandstop_8000_12000hz.txt\n";
+    std::cout << "- fir_lowpass_hann_10000hz.txt\n";
+    std::cout << "- fir_lowpass_hamming_10000hz.txt\n";
+    std::cout << "- fir_lowpass_blackman_10000hz.txt\n";
+    std::cout << "- plot_fir_responses.py\n\n";
+    std::cout << "Run: python3 plot_fir_responses.py (requires numpy, matplotlib, scipy)\n";
+}