Compactcodes.m

Inputfile = '';

% Chorus function (Effector)
function [outputSignal, Fs] = chorusOrig(inputFile)
if nargin < 1
    inputFile = 'inputfiles/solo_F.wav';
end

outputFile = 'outputs/chorusOrig.wav';

numVoices = 21;
pitchShifts = -8.0:0.8:8.0;  % semitones
delays = 0.1:-0.005:0;        % seconds

windowLength = 1024;
overlap = round(0.75 * windowLength);

[x, fs] = audioread(inputFile);
if size(x, 2) > 1
    x = mean(x, 2);
end
x = x(:);

y = zeros(size(x));

for v = 1:numVoices
    fprintf('Processing voice %d (shift: %.1f semitones)\n', v, pitchShifts(v));

    shifted = pitchShift(x, fs, pitchShifts(v), windowLength, overlap);

    delaySamples = round(delays(v) * fs);
    if delaySamples >= 0
        delayed = [zeros(delaySamples, 1); shifted(1:end-delaySamples)];
    else
        delayed = [shifted(-delaySamples+1:end); zeros(-delaySamples, 1)];
    end

    y = y + delayed / numVoices;
end

y = y / max(abs(y)) * 0.95;
audiowrite(outputFile, y, fs);

outputSignal = y;
Fs = fs;
end

% Filter function (Effector)
function [outputSignal, Fs] = filtering(inputFile)
if nargin < 1
    inputFile = 'inputfiles/music-sample.wav';  % fallback
end

cutoffLow = 300;   % Hz
cutoffHigh = 3000; % Hz
[audio, fs] = audioread(inputFile);
X = fft(audio);
f = (0:length(audio)-1) * (fs / length(audio));  % Frequency vector
% bandpass filter
filter = (f >= cutoffLow) & (f <= cutoffHigh);
X_filtered = X .* filter';
% Inverse FFT
audio_filtered = real(ifft(X_filtered));
% Normalize
audio_filtered = audio_filtered / max(abs(audio_filtered));
audiowrite('outputs/filtered.wav', audio_filtered, fs);

outputSignal = audioread('outputs/filtered.wav');   
[~, Fs] = audioread(inputFile);   
end

%Gender conversion function (Effector)
function [outputSignal, Fs] = Gender_style_conversion(inputFile, targetStyle)
if nargin < 1
    inputFile = 'inputfiles/voice-sample.wav';
end
if nargin < 2
    targetStyle = "feminine";
end

outputFile = 'outputs/gender_style_smooth.wav';

if ~exist('outputs', 'dir')
    mkdir('outputs');
end

[x, fs] = audioread(inputFile);

if size(x, 2) > 1
    x = mean(x, 2);
end

x = x(:);
x = x / max(abs(x) + eps);

fprintf('Input duration: %.2f seconds\n', length(x) / fs);
fprintf('Sampling rate: %d Hz\n', fs);


switch targetStyle

    case "feminine"
        pitchShiftSemitones = 2.0;
        formantScale = 1.06;

    case "masculine"
        medianF0 = estimateMedianF0(x, fs);
        targetF0 = 145;
        pitchShiftSemitones = 12 * log2(targetF0 / medianF0);
        pitchShiftSemitones = max(min(pitchShiftSemitones, -1.5), -6.0);
        formantScale = 2^(pitchShiftSemitones / 30);
        formantScale = max(min(formantScale, 0.97), 0.86);
        fprintf('Estimated median F0: %.2f Hz\n', medianF0);
        fprintf('Target F0: %.2f Hz\n', targetF0);

    otherwise
        error('targetStyle must be "feminine" or "masculine".');
end

y_pitch = smoothPitchShiftMultiStage(x, fs, pitchShiftSemitones);
fprintf('After pitch shift: max = %.6f, RMS = %.6f\n', max(abs(y_pitch)), rms(y_pitch));

y_formant = smoothFormantStyle(y_pitch, fs, formantScale);
fprintf('After formant style: max = %.6f, RMS = %.6f\n', max(abs(y_formant)), rms(y_formant));

y_final = smoothVoiceEQ(y_formant, fs, targetStyle);
y_final = shortFadeInOut(y_final, fs, 0.015);

if max(abs(y_final)) < 1e-6
    error('Output is almost silent. Processing failed.');
end

y_final = y_final / max(abs(y_final) + eps) * 0.95;

audiowrite(outputFile, y_final, fs);
fprintf('Output saved to: %s\n', outputFile);

outputSignal = y_final;
Fs = fs;

end

function y = smoothFormantStyle(x, fs, formantScale)
    x = x(:);
    N = length(x);
    windowLength = 2048;
    overlap = round(0.75 * windowLength);
    win = hamming(windowLength, 'periodic');
    [S, ~, ~] = stft(x, fs, 'Window', win, 'OverlapLength', overlap, ...
        'FFTLength', windowLength, 'Centered', false);
    mag = abs(S);
    phase = angle(S);
    numBins = size(S, 1);
    numFrames = size(S, 2);
    magNew = zeros(size(mag));
    binIndex = (1:numBins)';
    smoothBins = 55;
    for n = 1:numFrames
        currentMag = mag(:, n);
        logMag = log(currentMag + 1e-8);
        envelope = movmean(logMag, smoothBins);
        detail = logMag - envelope;
        sourceIndex = (binIndex - 1) / formantScale + 1;
        sourceIndex(sourceIndex < 1) = 1;
        sourceIndex(sourceIndex > numBins) = numBins;
        warpedEnvelope = interp1(binIndex, envelope, sourceIndex, 'linear');
        gainLog = min(max((warpedEnvelope + detail) - logMag, log(0.55)), log(1.6));
        magNew(:, n) = exp(logMag + gainLog);
    end
    Y = magNew .* exp(1i * phase);
    y = istft(Y, fs, 'Window', win, 'OverlapLength', overlap, ...
        'FFTLength', windowLength, 'Centered', false);
    y = matchLength(real(y), N);
end

function y = smoothVoiceEQ(x, fs, targetStyle)
    x = x(:);
    N = length(x);
    X = fft(x);
    f = (0:N-1)' * fs / N;
    fMirror = min(f, fs - f);
    gain = ones(N, 1);
    switch targetStyle
        case "feminine"
            gain(fMirror < 90) = gain(fMirror < 90) * 0.75;
            gain(fMirror >= 2200 & fMirror <= 6000) = gain(fMirror >= 2200 & fMirror <= 6000) * 1.08;
            gain(fMirror > 9000) = gain(fMirror > 9000) * 0.88;
        case "masculine"
            gain(fMirror >= 90  & fMirror <= 280)  = gain(fMirror >= 90  & fMirror <= 280)  * 1.30;
            gain(fMirror > 280  & fMirror <= 750)  = gain(fMirror > 280  & fMirror <= 750)  * 1.15;
            gain(fMirror >= 1800 & fMirror <= 3500) = gain(fMirror >= 1800 & fMirror <= 3500) * 0.90;
            gain(fMirror > 3500) = gain(fMirror > 3500) * 0.78;
            gain(fMirror < 55)  = gain(fMirror < 55)  * 0.65;
    end
    y = real(ifft(X .* gain));
end

function y = shortFadeInOut(x, fs, fadeSeconds)
    y = x(:);
    fadeLength = min(round(fadeSeconds * fs), floor(length(y) / 2));
    y(1:fadeLength) = y(1:fadeLength) .* linspace(0, 1, fadeLength)';
    y(end-fadeLength+1:end) = y(end-fadeLength+1:end) .* linspace(1, 0, fadeLength)';
end

function y = matchLength(y, targetLength)
    y = y(:);
    if length(y) > targetLength
        y = y(1:targetLength);
    elseif length(y) < targetLength
        y = [y; zeros(targetLength - length(y), 1)];
    end
end

function medianF0 = estimateMedianF0(x, fs)
    x = x(:);
    frameLength = round(0.040 * fs);
    hop = round(0.010 * fs);
    minLag = floor(fs / 450);
    maxLag = ceil(fs / 80);
    numFrames = floor((length(x) - frameLength) / hop) + 1;
    f0List = [];
    win = hamming(frameLength, 'periodic');
    for i = 1:numFrames
        startIdx = (i - 1) * hop + 1;
        frame = x(startIdx:startIdx + frameLength - 1) .* win;
        if rms(frame) < 0.01 * rms(x)
            continue;
        end
        r = xcorr(frame);
        r = r(frameLength:end);
        [peakVal, peakIdx] = max(r(minLag:maxLag));
        if peakVal > 0.25 * r(1)
            lag = peakIdx + minLag - 1;
            f0 = fs / lag;
            if f0 >= 80 && f0 <= 450
                f0List(end+1, 1) = f0; %#ok<AGROW>
            end
        end
    end
    if isempty(f0List)
        warning('Could not estimate F0 reliably. Using fallback value 220 Hz.');
        medianF0 = 220;
    else
        medianF0 = median(f0List);
    end
end

function y = smoothPitchShiftMultiStage(x, fs, totalSemitones)
    x = x(:);
    numStages = max(1, ceil(abs(totalSemitones) / 1.5));
    stepShift = totalSemitones / numStages;
    y = x;
    fprintf('Applying pitch shift in %d stages of %.2f semitones each\n', numStages, stepShift);
    for s = 1:numStages
        y = smoothPitchShiftPV(y, fs, stepShift);
        if max(abs(y)) < 1e-7
            warning('Stage %d produced nearly silent output. Falling back.', s);
            y = x;
            break;
        end
        y = removeTinyClicksSafe(y, fs);
        fprintf('  Stage %d/%d complete | max = %.6f | RMS = %.6f\n', s, numStages, max(abs(y)), rms(y));
    end
    y = matchLength(y, length(x));
end

function y = smoothPitchShiftPV(x, fs, semitones)
    x = x(:);
    originalLength = length(x);
    alpha = 2^(semitones / 12);
    stretched = phaseVocoderTimeStretch(x, fs, alpha);
    if max(abs(stretched)) < 1e-8
        warning('Time-stretch produced near silence. Returning original signal.');
        y = x;
        return;
    end
    y = interp1((1:length(stretched))', stretched, ...
        linspace(1, length(stretched), originalLength)', 'pchip', 0);
    y = matchLength(y(:), originalLength);
    if max(abs(y)) > 1e-8
        y = y / max(abs(y) + eps) * max(abs(x));
    else
        warning('Pitch shifted output is near silent. Returning original signal.');
        y = x;
    end
end

function y = phaseVocoderTimeStretch(x, ~, stretchFactor)
    x = x(:);
    stretchFactor = max(min(stretchFactor, 1.25), 0.75);
    N = length(x);
    windowLength = 2048;
    analysisHop = windowLength / 4;
    synthesisHop = round(analysisHop * stretchFactor);
    win = hann(windowLength, 'periodic');
    numFrames = floor((N - windowLength) / analysisHop) + 1;
    if numFrames < 2
        y = x;
        return;
    end
    outputLength = synthesisHop * (numFrames - 1) + windowLength;
    y = zeros(outputLength, 1);
    winSum = zeros(outputLength, 1);
    omega = 2 * pi * (0:windowLength-1)' / windowLength;
    previousPhase = zeros(windowLength, 1);
    synthesisPhase = zeros(windowLength, 1);
    for frameIndex = 1:numFrames
        inputStart = (frameIndex - 1) * analysisHop + 1;
        frame = x(inputStart:inputStart + windowLength - 1) .* win;
        X = fft(frame);
        mag = abs(X);
        phase = angle(X);
        if frameIndex == 1
            synthesisPhase = phase;
        else
            deltaPhase = phase - previousPhase - omega * analysisHop;
            deltaPhase = deltaPhase - 2*pi*round(deltaPhase / (2*pi));
            synthesisPhase = synthesisPhase + (omega + deltaPhase / analysisHop) * synthesisHop;
        end
        previousPhase = phase;
        outputStart = (frameIndex - 1) * synthesisHop + 1;
        outputRange = outputStart:outputStart + windowLength - 1;
        frameOut = real(ifft(mag .* exp(1i * synthesisPhase))) .* win;
        y(outputRange) = y(outputRange) + frameOut;
        winSum(outputRange) = winSum(outputRange) + win.^2;
    end
    valid = winSum > 1e-8;
    y(valid) = y(valid) ./ winSum(valid);
    y = real(y);
    if max(abs(y)) > 1e-8
        y = y / max(abs(y) + eps) * max(abs(x));
    end
end

function y = removeTinyClicksSafe(x, fs)
    y = x(:);
    fadeLength = min(round(0.008 * fs), floor(length(y) / 2));
    if fadeLength > 1
        y(1:fadeLength) = y(1:fadeLength) .* linspace(0, 1, fadeLength)';
        y(end-fadeLength+1:end) = y(end-fadeLength+1:end) .* linspace(1, 0, fadeLength)';
    end
    smoothing = 0.04;
    for n = 2:length(y)
        y(n) = (1 - smoothing) * y(n) + smoothing * y(n-1);
    end
end

%Equalizer function (Effector)
function [outputSignal, Fs] = graphicEqualizer(inputFile)
if nargin < 1
    inputFile = 'inputfiles/music-sample.wav'; 
end

[audio, fs] = audioread(inputFile);
audio = audio(:,1);
audio = audio / max(abs(audio));

windowSize = 1024;
hopSize = windowSize / 4;

lowGain  = 1.5;   % boost bass
midGain  = 1.0;   % keep middle
highGain = 0.6;   % reduce treble

numWindows = floor((length(audio) - windowSize)/hopSize) + 1;
output = zeros(size(audio));

for i = 1:numWindows
    
    startIdx = (i-1)*hopSize + 1;
    frame = audio(startIdx:startIdx+windowSize-1);
    win = hamming(windowSize);
    frameWindowed = frame .* win;
    F = fft(frameWindowed);
    N = length(F);
    gain = ones(N,1);
    lowEnd  = floor(N * 0.2);
    midEnd  = floor(N * 0.6);
    gain(1:lowEnd) = lowGain;
    gain(lowEnd+1:midEnd) = midGain;
    gain(midEnd+1:end) = highGain;
    F_eq = F .* gain;
    frameOut = real(ifft(F_eq)); 
    frameOut = frameOut .* win;
    output(startIdx:startIdx+windowSize-1) = ...
        output(startIdx:startIdx+windowSize-1) + frameOut;
end

output = output / max(abs(output));

audiowrite('outputs/equalized_voice.wav', output, fs);


outputSignal = audioread('outputs/equalized_voice.wav');   
[~, Fs] = audioread(inputFile);   
end

%Robotic (effecter)
function [outputSignal, Fs] = roboticdistortion(inputFile)
if nargin < 1
    inputFile = 'inputfiles/voice-sample.wav';  % fallback
end

[audio, fs] = audioread(inputFile); 
audio = audio(:,1);  
audio = audio / max(abs(audio));  % Normalize
% Parameters
windowSize = 1024;      % FFT size 
hopSize = windowSize / 4;  % Overlap, 75%
shiftAmount = 0.001; %  how much we shift the frequencies, smaller=more natural
phaseRand = 0.3; % phase randomization amount
%noiseGateThreshold = 0.15; % decides how much noise to remove, higher = more aggressive
% Process with STFT 
numWindows = floor((length(audio) - windowSize)/hopSize) + 1;
output = zeros(size(audio));

for i = 1:numWindows
    startIdx = (i-1)*hopSize + 1;
    frame = audio(startIdx:startIdx+windowSize-1);
    % Apply window
    win = hamming(windowSize);
    frameWindowed = frame .* win;
    % DFT
    F = fft(frameWindowed);
    
    % ROBOT EFFECT 
    % Frequency shift, makes it sound mechanical
    shift = round(length(F) * shiftAmount);  % shift by shiftAmount% of spectrum
    F_robot = [F(shift+1:end); zeros(shift,1)];  % simple shift
    % Randomize phase for more metallic sound
    mag = abs(F_robot);
    % uncomment one of the two lines below for random or zero phase
    phase = angle(F_robot) + (rand(size(F_robot))-0.5)*phaseRand; % random phase
    %phase = zeros(size(F_robot));   % zero phase instead of random, for more robotic sound
    F_robot = mag .* exp(1j * phase);
    
    % Inverse DFT
    frameOut = real(ifft(F_robot));
    frameOut = frameOut .* win;  % window again
    % Overlap-add
    output(startIdx:startIdx+windowSize-1) = output(startIdx:startIdx+windowSize-1) + frameOut;
end
% Normalize
output = output / max(abs(output));
% Save output
audiowrite('outputs/robotic_voice.wav', output, fs);


outputSignal = audioread('outputs/robotic_voice.wav');   
[~, Fs] = audioread(inputFile);   
end

%Paramter sensitivity test
% DFT Size Comparison — bar spectra at 10 Hz resolution
clear; clc; close all;

%% --- PARAMETERS ---
filename   = 'inputfiles/music-sample.wav';
window_end = 2.0;    % seconds to analyse
max_freq   = 4000;   % Hz to display on x-axis
bin_width  = 10;     % Hz per bar
hop        = 256;    % hop size (fixed across all DFT sizes)

%% --- LOAD & CROP ---
[x, fs] = audioread(filename);
x = mean(x, 2);
x = x(1 : min(end, round(window_end * fs)));

%% --- STFT RECONSTRUCT FOR EACH DFT SIZE ---
function x_out = stft_reconstruct(x, N, hop)
    num_frames = floor((length(x) - N) / hop) + 1;
    x_out      = zeros(length(x), 1);
    win        = hann(N);
    for f = 1:num_frames
        i1    = (f-1)*hop + 1;
        i2    = i1 + N - 1;
        frame = x(i1:i2) .* win;
        X     = fft(frame, N);
        rec   = real(ifft(X, N)) .* win;   % overlap-add
        x_out(i1:i2) = x_out(i1:i2) + rec;
    end
    x_out = x_out / max(abs(x_out) + 1e-10);
end

%% --- COMPUTE AVERAGE MAGNITUDE PER 10 Hz BIN ---
function mag_bins = compute_bins(x, fs, N, bin_width, max_freq)
    X        = abs(fft(x, N));
    X        = X(1:floor(N/2));
    freqs    = (0:floor(N/2)-1) * fs/N;
    edges    = 0 : bin_width : max_freq;
    mag_bins = zeros(1, length(edges)-1);
    for i = 1:length(edges)-1
        idx = freqs >= edges(i) & freqs < edges(i+1);
        if any(idx), mag_bins(i) = mean(X(idx)); end
    end
end

%% --- PROCESS & SAVE ---
sizes = {4096, 2048, 1024, 512};
recons = cell(1,4);
for k = 1:4
    N          = sizes{k};
    recons{k}  = stft_reconstruct(x, N, hop);
    audiowrite(sprintf('outputs/output_N%d.wav', N), recons{k}, fs);
end

%% --- COMPUTE BINS ---
centers = (bin_width/2 : bin_width : max_freq - bin_width/2);
shift   = 5;
bins = cell(1,4);
for k = 1:4
    bins{k} = compute_bins(recons{k}, fs, sizes{k}, bin_width, max_freq);
end

%% --- PLOT ---
labels = {'4096','2048','1024','512'};
colors = [0.2 0.5 0.9;
          0.9 0.4 0.2;
          0.2 0.8 0.4;
          0.8 0.2 0.6];

figure('Position', [100 100 1100 750]);
for p = 1:3
    subplot(3,1,p);
    bar(centers,       bins{p},   0.8, 'FaceColor', colors(p,:),   'FaceAlpha', 0.5, 'EdgeColor', 'none'); hold on;
    bar(centers+shift, bins{p+1}, 0.8, 'FaceColor', colors(p+1,:), 'FaceAlpha', 0.5, 'EdgeColor', 'none');
    xlabel('Frequency (Hz)'); ylabel('Magnitude');
    title(sprintf('N = %s  vs  N = %s', labels{p}, labels{p+1}));
    legend(sprintf('N=%s', labels{p}), sprintf('N=%s (+%dHz)', labels{p+1}, shift), 'Location','northeast');
    xlim([0 max_freq]); grid on;
end
sgtitle('DFT Size Comparison — STFT reconstruct, avg magnitude per 10 Hz bin');

%% --- FIGURE 2: HOP SIZE COMPARISON (fixed N=16384) ---
N_fixed    = 8192;
hops       = {256, 512, 1024, 2048};
hop_labels = {'256','512','1024','2048'};
hop_colors = [0.3 0.7 0.9;
              0.9 0.6 0.1;
              0.4 0.85 0.5;
              0.75 0.3 0.75];

hop_recons = cell(1,4);
hop_bins   = cell(1,4);
for k = 1:4
    hop_recons{k} = stft_reconstruct(x, N_fixed, hops{k});
    audiowrite(sprintf('outputs/output_hop%d.wav', hops{k}), hop_recons{k}, fs);
    hop_bins{k}   = compute_bins(hop_recons{k}, fs, N_fixed, bin_width, max_freq);
end

figure('Position', [200 100 1100 750]);
for p = 1:3
    subplot(3,1,p);
    bar(centers,       hop_bins{p},   0.8, 'FaceColor', hop_colors(p,:),   'FaceAlpha', 0.5, 'EdgeColor', 'none'); hold on;
    bar(centers+shift, hop_bins{p+1}, 0.8, 'FaceColor', hop_colors(p+1,:), 'FaceAlpha', 0.5, 'EdgeColor', 'none');
    xlabel('Frequency (Hz)'); ylabel('Magnitude');
    title(sprintf('hop = %s  vs  hop = %s', hop_labels{p}, hop_labels{p+1}));
    legend(sprintf('hop=%s', hop_labels{p}), sprintf('hop=%s (+%dHz)', hop_labels{p+1}, shift), 'Location','northeast');
    xlim([0 max_freq]); grid on;
end
sgtitle('Hop Size Comparison — N=16384 fixed, avg magnitude per 10 Hz bin');

%% --- FIGURE 3: OBSERVATION WINDOW LENGTH COMPARISON ---
fs_fig3    = 22500;
win_lens   = {0.1, 0.5, 1.5, 3.0};
win_labels = {'0.1s','0.5s','1.5s','3.0s'};
win_colors = [0.2 0.6 0.9;
              0.9 0.3 0.3;
              0.3 0.8 0.4;
              0.85 0.6 0.1];

% resample audio to fs_fig3 and loop to fill if needed
x_rs = resample(x, fs_fig3, fs);
n_needed = round(20.0 * fs_fig3);
while length(x_rs) < n_needed
    x_rs = [x_rs; x_rs];
end
x_rs = x_rs(1:n_needed);

win_bins = cell(1,4);
for k = 1:4
    % crop to exact window length — this IS the full observation window
    N_w    = round(win_lens{k} * fs_fig3);
    x_crop = x_rs(1:N_w);
    df     = fs_fig3 / N_w;
    fprintf('Window=%.1fs  N_w=%d  Δf=%.4f Hz\n', win_lens{k}, N_w, df);
    audiowrite(sprintf('outputs/output_win%s.wav', win_labels{k}), x_crop, fs_fig3);
    win_bins{k} = compute_bins(x_crop, fs_fig3, N_w, bin_width, max_freq);
end

figure('Position', [300 100 1100 750]);
for p = 1:3
    subplot(3,1,p);
    N_lo  = round(win_lens{p}   * fs_fig3);
    N_hi  = round(win_lens{p+1} * fs_fig3);
    df_lo = fs_fig3 / N_lo;
    df_hi = fs_fig3 / N_hi;

    bar(centers,       win_bins{p},   0.8, 'FaceColor', win_colors(p,:),   'FaceAlpha', 0.5, 'EdgeColor', 'none'); hold on;
    bar(centers+shift, win_bins{p+1}, 0.8, 'FaceColor', win_colors(p+1,:), 'FaceAlpha', 0.5, 'EdgeColor', 'none');
    xlabel('Frequency (Hz)'); ylabel('Magnitude');
    title(sprintf('win=%s (Δf=%.3fHz)  vs  win=%s (Δf=%.3fHz)', ...
        win_labels{p}, df_lo, win_labels{p+1}, df_hi));
    legend(win_labels{p}, sprintf('%s (+%dHz)', win_labels{p+1}, shift), 'Location','northeast');
    xlim([0 max_freq]); grid on;
end
sgtitle('Observation Window Comparison — fs=22500Hz, full window, avg magnitude per 10 Hz bin');
%Noise robustness test
function noiseRobustnessAnalysis()
    outputDir = fullfile('results', 'noise_robustness');
    if ~exist(outputDir, 'dir')
        mkdir(outputDir);
    end

    snrLevels = [30, 20, 10, 0];  % dB

    tests = { ...
        struct('name','Robotic',  'func',@roboticdistortion, 'input','inputfiles/voice-sample.wav'), ...
        struct('name','Equalizer','func',@graphicEqualizer,  'input','inputfiles/music-sample.wav'), ...
        struct('name','Filter',   'func',@filtering,         'input','inputfiles/music-sample.wav'), ...
        struct('name','Chorus',   'func',@chorusOrig,                                          'input','inputfiles/singingm4a-sample.wav'), ...
        struct('name','GenderF',  'func',@(f) Gender_style_conversion(f, "feminine"),          'input','inputfiles/voice-sample.wav'), ...
        struct('name','GenderM',  'func',@(f) Gender_style_conversion(f, "masculine"),         'input','inputfiles/vocals_F.wav') ...
    };

    for t = 1:length(tests)
        test = tests{t};
        fprintf(' %s \n', test.name);

        [cleanSig, Fs] = audioread(test.input);
        cleanSig = mean(cleanSig, 2);
        cleanSig = cleanSig / max(abs(cleanSig) + eps);

        % Baseline
        cleanPath = fullfile(tempdir, sprintf('clean_%s.wav', test.name));
        audiowrite(cleanPath, cleanSig, Fs);
        cleanProc = runEffectQuiet(test.func, cleanPath);
        cleanProc = cleanProc(:) / max(abs(cleanProc) + eps);

        N = 8192;
        freq = (0:N/2-1) * (Fs / N);
        half = 1:N/2;

        Y_cleanOut = abs(fft(cleanProc, N));

        fig = figure('Position', [100 100 1000 600], 'Visible', 'off');
        ax = axes('Parent', fig);
        plot(ax, freq, 20*log10(Y_cleanOut(half)/max(Y_cleanOut)+eps), 'k', 'LineWidth', 1.5);
        hold(ax, 'on');
        legendEntries = {'Clean'};

        for i = 1:length(snrLevels)
            snrIn = snrLevels(i);
            noisySig = addNoise(cleanSig, snrIn);
            noisySig = noisySig / max(abs(noisySig) + eps);

            noisyPath = fullfile(tempdir, sprintf('noisy_%s_%02d.wav', test.name, snrIn));
            audiowrite(noisyPath, noisySig, Fs);
            noisyProc = runEffectQuiet(test.func, noisyPath);
            noisyProc = noisyProc(:) / max(abs(noisyProc) + eps);

            L = min(length(cleanProc), length(noisyProc));
            cp = cleanProc(1:L);
            np = noisyProc(1:L);

            Y_noisyOut = abs(fft(np, N));
            outSNR = 10*log10(sum(cp.^2) / (sum((cp - np).^2) + eps));
            fprintf('  Input SNR=%2d dB | Output SNR=%6.2f dB\n', snrIn, outSNR);

            plot(ax, freq, 20*log10(Y_noisyOut(half)/max(Y_noisyOut)+eps), 'LineWidth', 1);
            legendEntries{end+1} = sprintf('SNR_{in}=%d dB', snrIn); %#ok<AGROW>
        end

        title(ax, sprintf('%s - output magnitute with regard to frequency domain', test.name));
        xlabel(ax, 'Frequency (Hz)'); ylabel(ax, 'Magnitude (dB)');
        legend(ax, legendEntries); grid(ax, 'on'); xlim(ax, [0 Fs/2]); ylim(ax, [-80 5]);
        saveas(fig, fullfile(outputDir, sprintf('%s.png', test.name)));
        close(fig);
    end

    fprintf('\nDone. Figures saved in: %s\n', outputDir);
end

function noisy = addNoise(sig, snrDb)
    Psig = mean(sig.^2);
    Pn = Psig / 10^(snrDb/10);
    noisy = sig + sqrt(Pn) * randn(size(sig));
end

function out = runEffectQuiet(fn, inputPath)
    [~, out] = evalc('fn(inputPath)');
    try
        clear playsnd
    catch
    end
end

%Sampling and resolution
function samplingResolutionAnalysis()

    outputDir = fullfile('results', 'sampling_resolution');
    if ~exist(outputDir, 'dir')
        mkdir(outputDir);
    end

    % Integer downsampling factors relative to each file's native Fs.
    % factor 1 = full rate (the original / baseline reference).
    % Each extra factor = one more full (expensive) effect run per test, so
    % we keep just the original plus one representative heavy downsample.
    factors = [1, 4];

    tests = { ...
        struct('name','Robotic',  'func',@roboticdistortion, 'input','inputfiles/voice-sample.wav'), ...
        struct('name','Equalizer','func',@graphicEqualizer,  'input','inputfiles/music-sample.wav'), ...
        struct('name','Filter',   'func',@filtering,         'input','inputfiles/music-sample.wav'), ...
        struct('name','Chorus',   'func',@chorusOrig,                                          'input','inputfiles/singingm4a-sample.wav'), ...
        struct('name','GenderF',  'func',@(f) Gender_style_conversion(f, "feminine"),          'input','inputfiles/voice-sample.wav'), ...
        struct('name','GenderM',  'func',@(f) Gender_style_conversion(f, "masculine"),         'input','inputfiles/vocals_F.wav') ...
    };

    N = 8192;            % DFT size used for every spectrum
    half = 1:N/2;

    for t = 1:length(tests)
        test = tests{t};
        fprintf('\n=== %s ===\n', test.name);

        [cleanSig, Fs] = audioread(test.input);
        cleanSig = mean(cleanSig, 2);                       % force mono
        cleanSig = cleanSig / max(abs(cleanSig) + eps);

        procRef  = processAtRate(test.func, cleanSig, Fs);
        procRef  = procRef(:) / max(abs(procRef) + eps);
        Yref     = avgSpectrum(procRef, N);     % averaged over whole signal
        refDb    = 20*log10(Yref(half)/max(Yref) + eps);
        freqRef  = (0:N/2-1) * (Fs / N);

        fig = figure('Position', [100 100 1100 800], 'Visible', 'off');

        ax1 = subplot(2,1,1, 'Parent', fig);
        hold(ax1, 'on');
        legend1 = {};

       
        dsFreq = cell(1, length(factors));
        dsDb   = cell(1, length(factors));

        for k = 1:length(factors)
            M = factors(k);
            FsDs = round(Fs / M);

            if M == 1
                freq = freqRef;
                yDb  = refDb;
            else
                sigDs = resample(cleanSig, 1, M);           % polyphase + LPF
                procDs = processAtRate(test.func, sigDs, FsDs);
                procDs = procDs(:) / max(abs(procDs) + eps);
                freq = (0:N/2-1) * (FsDs / N);
                Y    = avgSpectrum(procDs, N);
                yDb  = 20*log10(Y(half)/max(Y) + eps);
            end

            dsFreq{k} = freq;
            dsDb{k}   = yDb;

            plot(ax1, freq, yDb, 'LineWidth', 1.2);
            if M == 1
                legend1{end+1} = sprintf('M=1  ORIGINAL (Fs=%d Hz, Ny=%d Hz)', ...
                                         Fs, round(Fs/2)); %#ok<AGROW>
            else
                legend1{end+1} = sprintf('M=%d  (Fs=%d Hz, Ny=%d Hz)', ...
                                         M, FsDs, round(FsDs/2)); %#ok<AGROW>
            end
        end

        title(ax1, sprintf(['%s  -  original (M=1) vs downsampled ' ...
            '(band ends at each Nyquist)'], test.name));
        xlabel(ax1, 'Frequency (Hz)'); ylabel(ax1, 'Magnitude (dB)');
        legend(ax1, legend1, 'Location', 'southwest');
        grid(ax1, 'on'); xlim(ax1, [0 Fs/2]); ylim(ax1, [-80 5]);

        % ============ Subplot 2: deviation from the original ==============
        % For each downsampled spectrum, interpolate the ORIGINAL onto the
        % same (finer) frequency grid over the shared band [0, Fs_ds/2] and
        % plot (downsampled - original) in dB. Flat ~0 means the spectral
        % conclusions are unchanged by downsampling.
        ax2 = subplot(2,1,2, 'Parent', fig);
        hold(ax2, 'on');
        legend2 = {};

        for k = 2:length(factors)                            % skip M=1 (==0)
            freq = dsFreq{k}(:);
            refOnGrid = interp1(freqRef, refDb, freq, 'linear');
            diffDb = dsDb{k}(:) - refOnGrid(:);              % force column: avoid N/2 x N/2 broadcast
            plot(ax2, freq, diffDb, 'LineWidth', 1.0);
            FsDs = round(Fs / factors(k));
            legend2{end+1} = sprintf('M=%d  (vs original, band to %d Hz)', ...
                                     factors(k), round(FsDs/2)); %#ok<AGROW>
        end
        yline(ax2, 0, 'k--');

        title(ax2, sprintf(['%s  -  deviation of downsampled spectrum ' ...
            'from the original (dB)'], test.name));
        xlabel(ax2, 'Frequency (Hz)'); ylabel(ax2, '\Delta Magnitude (dB)');
        legend(ax2, legend2, 'Location', 'southwest');
        grid(ax2, 'on'); xlim(ax2, [0 Fs/2]); ylim(ax2, [-40 40]);

        saveas(fig, fullfile(outputDir, sprintf('%s.png', test.name)));
        close(fig);
    end

    fprintf('\nDone. Figures saved in: %s\n', outputDir);
end

% --- run an effect on an in-memory signal at a given sample rate ----------
function out = processAtRate(fn, sig, Fs)
    p = fullfile(tempdir, sprintf('sr_%d_%d.wav', Fs, round(rand*1e6)));
    audiowrite(p, sig / max(abs(sig) + eps), Fs);
    [~, out] = evalc('fn(p)');
    try
        clear playsnd
    catch
    end
    if exist(p, 'file'); delete(p); end
end

% --- Welch-style averaged magnitude spectrum over the WHOLE signal --------
% Avoids the trap of fft(sig, N) using only the first N samples, which can
% be silent (vocoder latency / leading delays) and yield max(Y)=0 -> NaN.
function mag = avgSpectrum(sig, N)
    sig = sig(:);
    if numel(sig) < N
        sig = [sig; zeros(N - numel(sig), 1)];
    end
    win  = hann(N, 'periodic');
    hop  = N / 2;                                   % 50% overlap
    nFr  = 1 + floor((numel(sig) - N) / hop);
    acc  = zeros(N, 1);
    for i = 1:nFr
        s     = (i-1)*hop + 1;
        frame = sig(s:s+N-1) .* win;
        acc   = acc + abs(fft(frame));
    end
    mag = acc / max(nFr, 1);
end