#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <algorithm>
#include <random> // For random number generation
#include <cstdint> // For fixed-size integer types

const int NUM_FRAMES = 256;
const int SAMPLES_PER_FRAME = 2048;
const float SAMPLE_RATE = 44100.0f;
const double TWO_PI = 6.28318530718;

// --- LFO Types ---
enum class LFOType {
    SINE,
    COSINE,
    SAW_UP,
    SAW_DOWN,
    TRIANGLE, // Added triangle wave
    SQUARE,    // Added square wave
    NUM_LFO_TYPES // Keep this at the end to easily get the count
};

// --- LFO Function ---
float lfo(float phase, LFOType type) {
    switch (type) {
        case LFOType::SINE:
            return std::sin(phase);
        case LFOType::COSINE:
            return std::cos(phase);
        case LFOType::SAW_UP:
            return fmod(phase, TWO_PI) / TWO_PI;
        case LFOType::SAW_DOWN:
            return 1.0f - (fmod(phase, TWO_PI) / TWO_PI);
        case LFOType::TRIANGLE: {
            float value = (2.0f * fmod(phase, TWO_PI) / TWO_PI) - 1.0f;
            return (value > 0) ? 1.0f - value : value + 1.0f;  // Convert to -1 to 1 range
        }
        case LFOType::SQUARE:
            return (fmod(phase, TWO_PI) < M_PI) ? 1.0f : -1.0f; // Simple square wave

        default:
            return 0.0f; // Should never happen, but good practice
    }
}



// --- FM Operator Function ---
float fmOperator(float time, float frequency, float modulationIndex, float modulator, float phaseOffset = 0.0f) {
    return std::sin(TWO_PI * frequency * time + modulationIndex * modulator + phaseOffset);
}

// --- Wavetable Generation Function ---
std::vector<std::vector<float>> generateFMWavetable(
    float baseFrequency,
    const std::vector<std::pair<float, float>>& operators,  // Frequency, Modulation Index
    const std::vector<std::tuple<LFOType, float, float>>& lfos = {}  // type, frequency, amplitude
) {
    std::vector<std::vector<float>> wavetable(NUM_FRAMES, std::vector<float>(SAMPLES_PER_FRAME, 0.0f));

    for (int frame = 0; frame < NUM_FRAMES; ++frame) {
        for (int sample = 0; sample < SAMPLES_PER_FRAME; ++sample) {
            float time = static_cast<float>(sample) / SAMPLES_PER_FRAME;
            float frameProgress = static_cast<float>(frame) / (NUM_FRAMES - 1); // 0 to 1

            // Apply LFOs
            float lfoValue = 0.0;
             for (const auto& lfo_params : lfos)
             {
                LFOType type;
                float lfoFreq, lfoAmp;
                std::tie(type, lfoFreq, lfoAmp) = lfo_params;
                lfoValue += lfoAmp * lfo(TWO_PI * lfoFreq * frameProgress, type);  //LFO applied to frames
             }

            // Start with the base frequency
            float currentSample = 0.0f;

            // Chain the operators for FM synthesis
            float modulator = 0.0; // Initialize modulator
            for (size_t i = 0; i < operators.size(); ++i)
            {
                float opFreq = operators[i].first;
                float opModIndex = operators[i].second;

                //apply lfoValue to either frequency or index
                opFreq = opFreq * (1 + lfoValue); // Example frequency modulation by LFO
                //opModIndex = opModIndex * (1.0 + lfoValue); // or modulate modulation index
                if (i == 0) //first operator is modulated by time.
                {
                   modulator = fmOperator(time, opFreq, opModIndex, 0); //first op can only have phase mod
                }
                else
                {
                   //apply previous modulator value
                   modulator = fmOperator(time, opFreq, opModIndex, modulator);
                }

                if (i == operators.size() - 1) //if we've reach the final op, assign the output value
                {
                   currentSample = modulator;
                }
            }

            wavetable[frame][sample] = currentSample;
        }
    }
    return wavetable;
}
// --- Wavetable Generation Function (with Controlled Randomization) ---
std::vector<std::vector<float>> generateRandomFMWavetable() {
    std::vector<std::vector<float>> wavetable(NUM_FRAMES, std::vector<float>(SAMPLES_PER_FRAME, 0.0f));

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> base_freq_dist(50.0, 400.0);  // Narrower base frequency range
    std::uniform_real_distribution<> freq_ratio_dist(0.5, 4.0);     // Frequency ratio (around harmonic series)
    std::uniform_real_distribution<> mod_index_dist(0.0, 10.0);    // Reduced max modulation index
    std::uniform_real_distribution<> lfo_freq_dist(0.005, 0.5);    // Slower LFOs for frame modulation
    std::uniform_real_distribution<> lfo_amp_dist(0.0, 0.5);      // Reduced LFO amplitude
    std::uniform_int_distribution<> lfo_type_dist(0, static_cast<int>(LFOType::NUM_LFO_TYPES) - 1);
    std::uniform_int_distribution<> num_operators_dist(2, 4);     // Slightly fewer operators
    std::uniform_int_distribution<> num_lfos_dist(0, 2);          // Fewer LFOs

    float baseFrequency = base_freq_dist(gen);
    int numOperators = num_operators_dist(gen);
    int numLFOs = num_lfos_dist(gen);

    std::vector<std::pair<float, float>> operators;
    for (int i = 0; i < numOperators; ++i) {
        // Key Change: Operator frequencies are *ratios* of the base frequency
        float freqRatio = freq_ratio_dist(gen);
        // Bias towards harmonic ratios (integers and simple fractions)
        if (gen() % 2 == 0) { // 50% chance of being near a harmonic
            freqRatio = std::round(freqRatio * 2.0) / 2.0; // Round to nearest 0.5
        }
        operators.emplace_back(baseFrequency * freqRatio, mod_index_dist(gen));
    }

    std::vector<std::tuple<LFOType, float, float>> lfos;
    for (int i = 0; i < numLFOs; ++i) {
        lfos.emplace_back(
            static_cast<LFOType>(lfo_type_dist(gen)),
            lfo_freq_dist(gen),
            lfo_amp_dist(gen)
        );
    }

    for (int frame = 0; frame < NUM_FRAMES; ++frame) {
        for (int sample = 0; sample < SAMPLES_PER_FRAME; ++sample) {
            float time = static_cast<float>(sample) / SAMPLES_PER_FRAME;
            float frameProgress = static_cast<float>(frame) / (NUM_FRAMES - 1);

            float lfoValue = 0.0;
            for (const auto& lfo_params : lfos) {
                LFOType type;
                float lfoFreq, lfoAmp;
                std::tie(type, lfoFreq, lfoAmp) = lfo_params;
                lfoValue += lfoAmp * lfo(TWO_PI * lfoFreq * frameProgress, type);
            }

            float currentSample = 0.0f;
            float modulator = 0.0;
            for (size_t i = 0; i < operators.size(); ++i) {
                float opFreq = operators[i].first;
                float opModIndex = operators[i].second;
                opFreq = opFreq * (1.0 + lfoValue);  // Apply LFO to frequency
                if (i == 0) {
                    modulator = fmOperator(time, opFreq, opModIndex, 0);
                } else {
                    modulator = fmOperator(time, opFreq, opModIndex, modulator);
                }
                if (i == operators.size() - 1) {
                    currentSample = modulator;
                }
            }
            wavetable[frame][sample] = currentSample;
        }
    }
    return wavetable;
}
void writeWavetableToFile(const std::vector<std::vector<float>>& wavetable, const char* filename) {
    std::ofstream outputFile(filename, std::ios::binary);
    if (!outputFile) {
        std::cerr << "Error opening file for writing." << std::endl;
        return;
    }

    // Corrected hex header as an array of unsigned chars
    unsigned char header[] = {
        0x52, 0x49, 0x46, 0x46, 0x80, 0x00, 0x20, 0x00, 0x57, 0x41, 0x56, 0x45, 0x4A, 0x55, 0x4E, 0x4B,
        0x1C,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,
        0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x66,0x6D,0x74,0x20,0x10,0x00,
        0x00, 0x00, 0x03, 0x00, 0x01, 0x00,
        0x44, 0xAC, 0x00, 0x00, 0x10, 0xB1, 0x02, 0x00, 0x04, 0x00, 0x20, 0x00, 0x63, 0x6C, 0x6D, 0x20,
        0x30, 0x00, 0x00, 0x00, 0x3C, 0x21, 0x3E, 0x32, 0x30, 0x34, 0x38, 0x20, 0x30, 0x30, 0x30, 0x30,
        0x30, 0x30, 0x30, 0x30, 0x20, 0x77, 0x61, 0x76, 0x65, 0x74, 0x61, 0x62, 0x6C, 0x65, 0x20, 0x28,
        0x77, 0x77, 0x77, 0x2E, 0x78, 0x66, 0x65, 0x72, 0x72, 0x65, 0x63, 0x6F, 0x72, 0x64, 0x73, 0x2E,
        0x63, 0x6F, 0x6D, 0x29, 0x64, 0x61, 0x74, 0x61, 0x00, 0x00, 0x20, 0x00, 0x00, 0x00, 0x00, 0x00,
        0x00, 0x00, 0x00, 0x00
    };

    outputFile.write(reinterpret_cast<const char*>(header), sizeof(header));

    // Write audio data
    for (const auto& frame : wavetable) {
        outputFile.write(reinterpret_cast<const char*>(frame.data()), frame.size() * sizeof(float));
    }
    outputFile.close();
}

// --- Function to Check for Repetition ---
bool isWavetableRepeating(const std::vector<std::vector<float>>& wavetable, float tolerance = 1e-4f) {
    if (wavetable.empty() || wavetable[0].empty()) {
        return false; // Empty wavetable is not considered repeating
    }

    const auto& firstFrame = wavetable[0];
    for (size_t frame = 1; frame < wavetable.size(); ++frame) {
        float maxDiff = 0.0f;
        for (size_t sample = 0; sample < firstFrame.size(); ++sample) {
            maxDiff = std::max(maxDiff, std::abs(wavetable[frame][sample] - firstFrame[sample]));
        }
        if (maxDiff > tolerance) {
            return false; // Significant difference found, not repeating
        }
    }
    return true; // All frames are very similar to the first frame
}
// --- Function to Calculate Spectral Flatness ---
float calculateSpectralFlatness(const std::vector<std::vector<float>>& wavetable) {
    if (wavetable.empty() || wavetable[0].empty()) {
        return 0.0f; // Handle empty wavetable
    }

    std::vector<float> spectrum;
    for (const auto& frame : wavetable)
    {
        //perform fft on frame, and take the magnitudes.
        //in a real scenario, we'd use an FFT library like FFTW
        //here, we're going to *approximate* the spectrum using a
        //Discrete Cosine Transform (DCT-II), which is simpler to
        //implement and gives us a reasonable idea of frequency content
        //for real-valued signals.
        for (int k = 0; k < SAMPLES_PER_FRAME; ++k)
        {
            float sum = 0.0;
            for (int n = 0; n < SAMPLES_PER_FRAME; ++n)
            {
                sum += frame[n] * cos(M_PI / SAMPLES_PER_FRAME * (n + 0.5) * k);
            }
            spectrum.push_back(std::abs(sum)); // Magnitude
        }
    }

    // Ensure no zero values in the spectrum for geometric mean calculation
    for (float& val : spectrum) {
        if (val <= 0.0f) {
            val = 1e-6f;  // Replace with a small positive value
        }
    }

    // Calculate geometric mean
    double geometricMean = 1.0;
    for (float val : spectrum) {
        geometricMean *= val;
    }
    geometricMean = pow(geometricMean, 1.0 / spectrum.size());


    // Calculate arithmetic mean
    double arithmeticMean = std::accumulate(spectrum.begin(), spectrum.end(), 0.0) / spectrum.size();

    // Calculate spectral flatness (avoid division by zero)
    if (arithmeticMean <= 0.0) {
        return 0.0f;  // If arithmetic mean is zero or negative, flatness is zero
    }
    return static_cast<float>(geometricMean / arithmeticMean);
}

// --- Function to Calculate Average Absolute Difference (AAD) ---
float calculateAverageAbsoluteDifference(const std::vector<std::vector<float>>& wavetable) {
    if (wavetable.empty() || wavetable[0].empty()) {
        return 0.0f; // Handle empty wavetable
    }

    double totalDifference = 0.0;
    long long numDifferences = 0; // Use long long to prevent potential overflow

    for (const auto& frame : wavetable) {
        for (size_t i = 0; i < frame.size() - 1; ++i) {
            totalDifference += std::abs(frame[i + 1] - frame[i]);
            numDifferences++;
        }
    }
    // Avoid division by zero
    if (numDifferences == 0) {
        return 0.0f;
    }

    return static_cast<float>(totalDifference / numDifferences);
}

int main() {
    std::random_device rd;

    std::vector<std::vector<float>> randomWavetable;
    int maxAttempts = 100;
    int attempts = 0;
    float aadThreshold = 0.7f;

    do {
        randomWavetable = generateRandomFMWavetable();
        attempts++;
        if (attempts > maxAttempts) {
            std::cerr << "Failed to generate a suitable wavetable after " << maxAttempts << " attempts." << std::endl;
            return 1;
        }
        std::cout << "Attempt " <<std::to_string(attempts)<< "\n";
        // Check for repetition *and* excessive flatness (noise)
    } while (isWavetableRepeating(randomWavetable) || calculateAverageAbsoluteDifference(randomWavetable) > aadThreshold);

    writeWavetableToFile(randomWavetable, "random_wavetable.wav");
    std::cout << "Random, non-repeating, non-noisy wavetable created: random_wavetable.wav" << std::endl;
    std::cout << "Generated in " << attempts << " attempts." << std::endl;

    return 0;
}




int main_old() {
    // Define a carrier frequency for the examples.
    const float carrierFrequency = 440.0f;

   // Example 1: Simple 2-operator FM
    std::vector<std::pair<float, float>> operators1 = {
        {carrierFrequency, 0.0f},       // Carrier: freq, mod index
        {carrierFrequency * 0.5f, 5.0f}  // Modulator: freq, mod index
    };
    auto wavetable1 = generateFMWavetable(carrierFrequency, operators1);
    writeWavetableToFile(wavetable1, "wavetable_2op.wav");



    // Example 2: 3-operator FM with an LFO on the modulator frequency
    std::vector<std::pair<float, float>> operators2 = {
        {110.0f, 0.0f},       // Carrier (A2)
        {220.0f, 10.0f},     // Modulator 1
        {55.0f, 20.0f}      // Modulator 2
    };
    //LFO, freq, amp
    std::vector<std::tuple<LFOType, float, float>> lfos2 = {
        {LFOType::SINE, 0.5f, 0.2f} // Sine LFO at 0.5 Hz, amplitude 0.2
    };
    auto wavetable2 = generateFMWavetable(110.0f, operators2, lfos2); // Base freq of A2
    writeWavetableToFile(wavetable2, "wavetable_3op_lfo.wav");


    // Example 3:  More complex FM with multiple LFOs
    std::vector<std::pair<float, float>> operators3 = {
        {220.0f, 0.0f},   // Carrier
        {440.0f, 5.0f},   // Mod 1
        {110.0f, 10.0f},   // Mod 2
        {55.0f, 15.0f}    // Mod 3
    };

    std::vector<std::tuple<LFOType, float, float>> lfos3 = {
        {LFOType::SAW_UP, 1.0f, 0.3f},    // Saw up LFO (affects frame progression)
        {LFOType::TRIANGLE, 0.2f, 0.1f},   // Triangle LFO
        {LFOType::COSINE, 0.7f, 0.25f}     // Cosine LFO
    };

    auto wavetable3 = generateFMWavetable(220.0f, operators3, lfos3);
    writeWavetableToFile(wavetable3, "wavetable_complex_lfo.wav");

	// Example 4:  Simple sine wavetable
    std::vector<std::pair<float, float>> operators4 = {
        {440.0f, 0.0f},  // A4
    };
    auto wavetable4 = generateFMWavetable(440.0f, operators4);
    writeWavetableToFile(wavetable4, "wavetable_sine.wav");

	// Example 5: Test all LFO types.  Use a single operator, but modulate it with ALL the LFOs
	std::vector<std::pair<float, float>> operators5 = {
	   {110.0f, 0.0f} // A2
	};

	std::vector<std::tuple<LFOType, float, float>> lfos5 = {
		{LFOType::SINE, 1.0f, 0.2f},       // Sine
		{LFOType::COSINE, 0.5f, 0.3f},    // Cosine
		{LFOType::SAW_UP, 0.2f, 0.4f},   // Saw Up
		{LFOType::SAW_DOWN, 0.1f, 0.5f},  // Saw Down
		{LFOType::TRIANGLE, 0.8f, 0.15f},  // Triangle
		{LFOType::SQUARE, 2.0f, 0.25f}     // Square
	};

	auto wavetable5 = generateFMWavetable(110.0f, operators5, lfos5);
	writeWavetableToFile(wavetable5, "wavetable_all_lfos.wav");

    std::cout << "Wavetable files created successfully." << std::endl;

    return 0;
}