Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
675 lines (549 sloc) 36.2 KB
% Matchering DSP Module 0.86 by SOUND.TOOLS
% https://github.com/SOUNDTOOLS
%
% You need to place mrswatson64.exe and activated brickwall limiter VST DLL in the folder specified in pathBin variable!
%
% All non-44100 audio will be resampled to 44100 for now
% Only stereo and mono supported for now
%
% This module made in procedural style, not OOP, sorry :)
%
% ----------------------------------------------------------------
function currentSessionStatus = matchering(pathWorkspace, fileTarget, fileReference, pathTemp, pathBin, fileBrickwallLimiter, maxLength, timeArea, fftSize, linLogOversampling, peakCompensationSteps)
disp('Matchering DSP Module 0.86 by SOUND.TOOLS')
disp('----------------------------------------------------------------')
currentSessionStatus = 2;
disp(['Target file: ', fileTarget])
disp(['Reference file: ', fileReference])
currentSessionStatus = 3;
filesIsOk = true;
%Loading the target
disp(['Loading target ', pathWorkspace, fileTarget, '...'])
try
[target, targetSamplingFrequency] = audioread([pathWorkspace, fileTarget]);
disp('Target file is loaded')
catch
disp('ERROR: Unknown audio stream in target file')
filesIsOk = false;
currentSessionStatus = 31;
end
%Loading the reference
disp(['Loading reference ', pathWorkspace, fileReference, '...'])
try
[reference, referenceSamplingFrequency] = audioread([pathWorkspace, fileReference]);
disp('Reference file is loaded')
catch
if filesIsOk == true
disp('ERROR: Unknown audio stream in reference file')
filesIsOk = false;
currentSessionStatus = 32;
else
disp('ERROR: Unknown audio stream in reference file too')
currentSessionStatus = 30;
end
end
if filesIsOk == true
% Hard link to 44100!
if targetSamplingFrequency ~= 44100
disp('Target track is not 44100! Resampling...')
target = resample(target, 44100, targetSamplingFrequency);
targetSamplingFrequency = 44100;
end
if referenceSamplingFrequency ~= 44100
disp('Reference track is not 44100! Resampling...')
reference = resample(reference, 44100, referenceSamplingFrequency);
referenceSamplingFrequency = 44100;
end
% Checking length of the target
maxLengthIsOk = true;
minLengthIsOk = true;
targetSize = size(target, 1);
disp(['Size of target WAV: ', num2str(targetSize) ,' samples or ', num2str(targetSize / targetSamplingFrequency), ' seconds'])
if targetSize > (double(maxLength) * targetSamplingFrequency * 60)
disp(['ERROR: Target track length is more than maximum (', num2str(double(maxLength)), ' minutes)'])
maxLengthIsOk = false;
currentSessionStatus = 34;
end
if targetSize < double(fftSize)
disp(['ERROR: Target track length is less than minimum (', num2str(double(fftSize)), ' samples)'])
minLengthIsOk = false;
currentSessionStatus = 40;
end
% Checking length of the reference
referenceSize = size(reference, 1);
disp(['Size of reference WAV: ', num2str(referenceSize), ' samples or ',num2str(referenceSize/referenceSamplingFrequency), ' seconds'])
if referenceSize > (double(maxLength) * referenceSamplingFrequency * 60)
if minLengthIsOk == true
if maxLengthIsOk == true
disp(['ERROR: Reference track length is more than maximum (', num2str(double(maxLength)), ' minutes)'])
maxLengthIsOk = false;
currentSessionStatus = 35;
else
disp(['ERROR: Reference track length is more than maximum (', num2str(double(maxLength)), ' minutes) too'])
currentSessionStatus = 33;
end
else
disp(['ERROR: Reference track length is more than maximum (', num2str(double(maxLength)),' minutes)'])
maxLengthIsOk = false;
currentSessionStatus = 42;
end
end
if referenceSize < double(fftSize)
if maxLengthIsOk == true
if minLengthIsOk == true
disp(['ERROR: Reference track length is less than minimum (', num2str(double(fftSize)), ' samples)'])
minLengthIsOk = false;
currentSessionStatus = 41;
else
disp(['ERROR: Reference track length is less than minimum (', num2str(double(fftSize)), ' samples) too'])
currentSessionStatus = 39;
end
else
disp(['ERROR: Reference track length is less than minimum (', num2str(double(fftSize)), ' samples)'])
minLengthIsOk = false;
currentSessionStatus = 43;
end
end
if maxLengthIsOk == true
if minLengthIsOk == true
areaSize = double(timeArea) * targetSamplingFrequency;
disp(['Max size of analyzed part: ', num2str(areaSize), ' samples or ', num2str(double(timeArea)), ' seconds'])
% Converting mono to stereo
if size(target, 2) == 1
disp('Target audio is mono. Converting it to stereo...')
target(:, 2) = target(:, 1);
end
if size(reference, 2) == 1
disp('Reference audio is mono. Converting it to stereo...')
reference(:, 2) = reference(:, 1);
end
% Determining audio files with more than 2 channels
channelsIsOk = true;
if size(target, 2) > 2
disp('ERROR: Target track has more than 2 channels')
channelsIsOk = false;
currentSessionStatus = 37;
end
if size(reference, 2) > 2
if channelsIsOk == true
disp('ERROR: Reference track has more than 2 channels')
channelsIsOk = false;
currentSessionStatus = 38;
else
disp('ERROR: Reference track has more than 2 channels too')
currentSessionStatus = 36;
end
end
if channelsIsOk == true
% Checking peaks and limiters
targetSameMaximumPoints = sum(target(:) == max(max(abs(target)))) + sum(target(:) == -max(max(abs(target))));
referenceSameMaximumPoints = sum(reference(:) == max(max(abs(reference)))) + sum(reference(:) == -max(max(abs(reference))));
if targetSameMaximumPoints > 8
if max(max(abs(target))) == 1
disp('WARNING: Audio clipping detected in the target! It is very recommended to use version without clipping!')
else
if targetSameMaximumPoints > 128
disp('WARNING: Applied limiter detected in the target! It is very recommended to use version without limiter!')
end
end
end
if referenceSameMaximumPoints > 8
if max(max(abs(reference))) == 1
disp('INFO: Audio clipping detected in the reference')
end
end
disp('--------------------------------')
disp('Stage 0 - Normalizing Reference')
currentSessionStatus = 4;
% Stage 0 - Normalizing Reference
LIMITED_MAXIMUM_POINT = 0.998138427734375;
if max(max(abs(reference))) >= LIMITED_MAXIMUM_POINT
finalAmplitudeCoefficient = 1;
disp('Reference is not changed. No final amplitude coefficient')
else
finalAmplitudeCoefficient = max(max(abs(reference))) / LIMITED_MAXIMUM_POINT;
if finalAmplitudeCoefficient < 0.000001
finalAmplitudeCoefficient = 0.000001;
end
reference = reference ./ finalAmplitudeCoefficient;
disp(['Reference is normalized. Final aplitude coefficient for target audio will be: ', num2str(finalAmplitudeCoefficient)])
end
disp('Stage 0 completed')
disp('--------------------------------')
% Stage 1 - Matching RMS
disp('Stage 1 - Matching RMS')
disp('Calculating mid channel of target...');
targetMid(:) = (target(:, 1) + target(:, 2)) ./ 2;
disp('Calculating mid channel of reference...');
referenceMid(:) = (reference(:, 1) + reference(:, 2)) ./ 2;
disp('Calculating side channel of target...');
targetSide(:) = (target(:, 1) - target(:, 2)) ./ 2;
disp('Calculating side channel of reference...');
referenceSide(:) = (reference(:, 1) - reference(:, 2)) ./ 2;
disp('Unloading original reference from memory...')
clearvars reference;
targetDivisions = fix(targetSize / areaSize) + 1;
disp(['Target will be didived into ', num2str(targetDivisions), ' parts'])
referenceDivisions = fix(referenceSize / areaSize) + 1;
disp(['Reference will be didived into ', num2str(referenceDivisions), ' parts'])
% Finding sizes of parts
targetPartSize = fix(targetSize / targetDivisions);
disp(['One part of target is ', num2str(targetPartSize), ' samples or ', num2str(targetPartSize / targetSamplingFrequency), ' seconds'])
referencePartSize = fix(referenceSize / referenceDivisions);
disp(['One part of reference is ', num2str(referencePartSize), ' samples or ', num2str(referencePartSize / referenceSamplingFrequency), ' seconds'])
disp('Calculating RMSes of target parts...');
targetRmses = zeros(1, targetDivisions);
for i = 1 : 1 : targetDivisions
targetRmses(i) = rms(targetMid( (targetPartSize * (i - 1) + 1) : (targetPartSize * i) ));
end
disp('Calculating RMSes of reference parts...');
referenceRmses = zeros(1, referenceDivisions);
for j = 1 : 1 : referenceDivisions
referenceRmses(j) = rms(referenceMid( (referencePartSize * (j - 1) + 1) : (referencePartSize*j) ));
end
% Average RMSes of mid channels
targetAverageRms = rms(targetRmses);
referenceAverageRms = rms(referenceRmses);
disp(['Extracting parts from target audio with RMS more than average ', num2str(targetAverageRms), '...'])
% May add array preallocation here in future
j = 1;
for i = 1 : 1 : targetDivisions
if targetRmses(i) >= (targetAverageRms)
targetAnalyzedPartNumbers(j) = i;
j = j + 1;
end
end
targetMidAnalyzedPartsTemp = zeros(size(targetAnalyzedPartNumbers, 2), targetPartSize);
targetSideAnalyzedPartsTemp = zeros(size(targetAnalyzedPartNumbers, 2), targetPartSize);
for i = 1 : 1 : size(targetAnalyzedPartNumbers, 2)
targetMidAnalyzedPartsTemp(i, :) = targetMid( (targetPartSize * (targetAnalyzedPartNumbers(i) - 1) + 1) : (targetPartSize * targetAnalyzedPartNumbers(i)) );
targetSideAnalyzedPartsTemp(i, :) = targetSide( (targetPartSize * (targetAnalyzedPartNumbers(i) - 1) + 1) : (targetPartSize * targetAnalyzedPartNumbers(i)) );
end
targetAnalyzedRmses = zeros(1, size(targetAnalyzedPartNumbers, 2));
for i = 1 : 1 : size(targetAnalyzedPartNumbers, 2)
targetAnalyzedRmses(i) = targetRmses(targetAnalyzedPartNumbers(i));
end
disp(['Extracting parts from reference audio with RMS more than average ', num2str(referenceAverageRms), '...'])
% May add array preallocation here in future
j = 1;
for i = 1 : 1 : referenceDivisions
if referenceRmses(i) >= (referenceAverageRms)
referenceAnalyzedPartNumbers(j) = i;
j = j + 1;
end
end
referenceMidAnalyzedPartsTemp = zeros(size(referenceAnalyzedPartNumbers, 2), referencePartSize);
referenceSideAnalyzedPartsTemp = zeros(size(referenceAnalyzedPartNumbers, 2), referencePartSize);
for i = 1 : 1 : size(referenceAnalyzedPartNumbers, 2)
referenceMidAnalyzedPartsTemp(i, :) = referenceMid( (referencePartSize * (referenceAnalyzedPartNumbers(i) - 1) + 1) : (referencePartSize * referenceAnalyzedPartNumbers(i)) );
referenceSideAnalyzedPartsTemp(i, :) = referenceSide( (referencePartSize * (referenceAnalyzedPartNumbers(i) - 1) + 1) : (referencePartSize * referenceAnalyzedPartNumbers(i)) );
end
referenceAnalyzedRmses = zeros(1, size(referenceAnalyzedPartNumbers, 2));
for i = 1 : 1 : size(referenceAnalyzedPartNumbers, 2)
referenceAnalyzedRmses(i) = referenceRmses(referenceAnalyzedPartNumbers(i));
end
disp('Unloading mid channel of reference from memory...')
clearvars referenceMid;
disp('Unloading side channel of reference from memory...')
clearvars referenceSide;
targetMatchingRms = rms(targetAnalyzedRmses);
if targetMatchingRms < 0.000001
targetMatchingRms = 0.000001;
end
referenceMatchingRms = rms(referenceAnalyzedRmses);
rmsCoefficient = referenceMatchingRms / targetMatchingRms;
disp(['The RMS coefficient is: ', num2str(rmsCoefficient)])
disp('Modifying amplitudes of target...')
outputStage1(:, 1) = targetMid .* rmsCoefficient; % Left channel of this is mid
outputStage1(:, 2) = targetSide .* rmsCoefficient; % Right channel of this is side
disp('Modifying amplitudes of extracted target parts...')
targetMidAnalyzedParts = targetMidAnalyzedPartsTemp .* rmsCoefficient;
targetSideAnalyzedParts = targetSideAnalyzedPartsTemp .* rmsCoefficient;
disp('Unloading mid channel of target from memory...')
clearvars targetMid;
disp('Unloading side channel of target from memory...')
clearvars targetSide;
disp('Partly cleaning memory...')
clearvars -except currentSessionStatus outputStage1 finalAmplitudeCoefficient pathWorkspace fileTarget fileReference pathTemp pathBin fileBrickwallLimiter maxLength timeArea fftSize linLogOversampling peakCompensationSteps targetSamplingFrequency targetDivisions targetPartSize targetMidAnalyzedParts targetSideAnalyzedParts referencePartSize referenceMidAnalyzedPartsTemp referenceSideAnalyzedPartsTemp referenceMatchingRms;
disp('Stage 1 completed')
disp('--------------------------------')
% Stage 2 - Matching FR
disp('Stage 2 - Matching FR')
currentSessionStatus = 5;
targetMidFftCount = fix(targetPartSize / double(fftSize));
targetSideFftCount = targetMidFftCount;
referenceMidFftCount = fix(referencePartSize / double(fftSize));
referenceSideFftCount = referenceMidFftCount;
% ---MID---
disp(['Calculating ', num2str(targetMidFftCount * size(targetMidAnalyzedParts, 1)), ' FFTs of target mid...'])
targetMidAverageFft = zeros(size(targetMidAnalyzedParts, 1), double(fftSize));
for j = 1 : 1 : size(targetMidAnalyzedParts, 1)
targetMidFft = zeros(targetMidFftCount, double(fftSize));
for i = 1 : 1 : targetMidFftCount
targetMidFft(i, :) = abs(fft(targetMidAnalyzedParts(j, (double(fftSize) * (i - 1) + 1) : (double(fftSize) * i) )));
end
if size(targetMidFft, 1) == 1
targetMidAverageFft(j, :) = targetMidFft;
else
targetMidAverageFft(j, :) = sum(targetMidFft) ./ targetMidFftCount;
end
clearvars targetMidFft;
end
disp(['Calculating ', num2str(referenceMidFftCount * size(referenceMidAnalyzedPartsTemp, 1)), ' FFTs of reference mid...'])
referenceMidAverageFft = zeros(size(referenceMidAnalyzedPartsTemp, 1), double(fftSize));
for j = 1 : 1 : size(referenceMidAnalyzedPartsTemp, 1)
referenceMidFft = zeros(referenceMidFftCount, double(fftSize));
for i = 1 : 1 : referenceMidFftCount
referenceMidFft(i, :) = abs(fft(referenceMidAnalyzedPartsTemp(j, (double(fftSize) * (i - 1) + 1) : (double(fftSize) * i) )));
end
if size(referenceMidFft, 1) == 1
referenceMidAverageFft(j, :) = referenceMidFft;
else
referenceMidAverageFft(j, :) = sum(referenceMidFft) ./ referenceMidFftCount;
end
clearvars referenceMidFft;
end
disp('Calculating average FFT of target mid...')
if size(targetMidAverageFft, 1) == 1
targetMidCompletedFft = targetMidAverageFft;
else
targetMidCompletedFft = sum(targetMidAverageFft) ./ size(targetMidAnalyzedParts, 1);
end
disp('Calculating average FFT of reference mid...')
if size(referenceMidAverageFft, 1) == 1
referenceMidCompletedFft = referenceMidAverageFft;
else
referenceMidCompletedFft = sum(referenceMidAverageFft) ./ size(referenceMidAnalyzedPartsTemp, 1);
end
disp('Calculating mid matching FFT...')
targetMidCompletedFft(targetMidCompletedFft < 0.000001) = 0.000001;
midMatchingFft = referenceMidCompletedFft ./ targetMidCompletedFft;
disp('Partly cleaning memory...')
clearvars referenceMidFft* targetMidFft*;
disp('Filtering mid matching FFT...')
gridLinear = targetSamplingFrequency / 2 * linspace(0, 1, double(fftSize) / 2 + 1);
gridLogarithmic = targetSamplingFrequency / 2 * logspace(log10(4 / double(fftSize)), 0, (double(fftSize) / 2) * double(linLogOversampling) + 1); % Oversampling here
midMatchingFftLogarithmic = interp1(gridLinear, midMatchingFft(1 : double(fftSize) / 2 + 1), gridLogarithmic, 'spline');
midMatchingFftLogarithmicFiltered = smooth(midMatchingFftLogarithmic, 0.075,'loess')';
midMatchingFftFiltered = interp1(gridLogarithmic, midMatchingFftLogarithmicFiltered, gridLinear,'spline');
midMatchingFftFiltered(double(fftSize) / 2 + 2 : double(fftSize)) = midMatchingFftFiltered(double(fftSize) / 2 : -1 : 2);
midMatchingFftFiltered(1) = 0; % Remove DC offset, lol
midMatchingFftFiltered(2) = midMatchingFft(2);
midMatchingFftFiltered(double(fftSize)) = midMatchingFft(double(fftSize));
disp('Calculating mid FIR for matching EQ...')
midFirTmp = real(ifft(midMatchingFftFiltered));
midFir = ifftshift(midFirTmp);
midFirHann = midFir .* hann(double(fftSize))';
clearvars midFirTmp midFir;
% ---SIDE---
disp(['Calculating ', num2str(targetSideFftCount * size(targetSideAnalyzedParts, 1)), ' FFTs of target side...'])
targetSideAverageFft = zeros(size(targetSideAnalyzedParts, 1), double(fftSize));
for j = 1 : 1 : size(targetSideAnalyzedParts, 1)
targetSideFft = zeros(targetSideFftCount, double(fftSize));
for i = 1 : 1 : targetSideFftCount
targetSideFft(i, :) = abs(fft(targetSideAnalyzedParts(j, (double(fftSize) * (i - 1) + 1) : (double(fftSize) * i) )));
end
if size(targetSideFft, 1) == 1
targetSideAverageFft(j, :) = targetSideFft;
else
targetSideAverageFft(j, :) = sum(targetSideFft) ./ targetSideFftCount;
end
clearvars targetSideFft;
end
disp(['Calculating ', num2str(referenceSideFftCount * size(referenceSideAnalyzedPartsTemp, 1)), ' FFTs of reference side...'])
referenceSideAverageFft = zeros(size(referenceSideAnalyzedPartsTemp, 1), double(fftSize));
for j = 1 : 1 : size(referenceSideAnalyzedPartsTemp, 1)
referenceSideFft = zeros(referenceSideFftCount, double(fftSize));
for i = 1 : 1 : referenceSideFftCount
referenceSideFft(i, :) = abs(fft(referenceSideAnalyzedPartsTemp(j, (double(fftSize) * (i - 1) + 1) : (double(fftSize) * i) )));
end
if size(referenceSideFft, 1) == 1
referenceSideAverageFft(j, :) = referenceSideFft;
else
referenceSideAverageFft(j, :) = sum(referenceSideFft) ./ referenceSideFftCount;
end
clearvars referenceSideFft;
end
disp('Calculating average FFT of target side...')
if size(targetSideAverageFft, 1) == 1
targetSideCompletedFft = targetSideAverageFft;
else
targetSideCompletedFft = sum(targetSideAverageFft) ./ size(targetSideAnalyzedParts, 1);
end
disp('Calculating average FFT of reference side...')
if size(referenceSideAverageFft, 1) == 1
referenceSideCompletedFft = referenceSideAverageFft;
else
referenceSideCompletedFft = sum(referenceSideAverageFft) ./ size(referenceSideAnalyzedPartsTemp, 1);
end
disp('Calculating side matching FFT...')
targetSideCompletedFft(targetSideCompletedFft < 0.000001) = 0.000001;
sideMatchingFft = referenceSideCompletedFft ./ targetSideCompletedFft;
disp('Partly cleaning memory...')
clearvars referenceSideFft* targetSideFft*;
disp('Filtering side matching FFT...')
sideMatchingFftLogarithmic = interp1(gridLinear, sideMatchingFft(1 : double(fftSize) / 2 + 1), gridLogarithmic, 'spline');
sideMatchingFftLogarithmicFiltered = smooth(sideMatchingFftLogarithmic, 0.075, 'loess')';
sideMatchingFftFiltered = interp1(gridLogarithmic, sideMatchingFftLogarithmicFiltered, gridLinear, 'spline');
sideMatchingFftFiltered(double(fftSize) / 2 + 2 : double(fftSize)) = sideMatchingFftFiltered(double(fftSize) / 2 : -1 : 2);
sideMatchingFftFiltered(1) = 0; % Remove DC offset, lol
sideMatchingFftFiltered(2) = sideMatchingFft(2);
sideMatchingFftFiltered(double(fftSize)) = sideMatchingFft(double(fftSize));
disp('Calculating side FIR for matching EQ...')
sideFirTemp = real(ifft(sideMatchingFftFiltered));
sideFir = ifftshift(sideFirTemp);
sideFirHann = sideFir.*hann(double(fftSize))';
clearvars sideFirTemp sideFir;
disp('Convolving target with calculated FIRs...')
tic;
outputStage1MidTemp = outputStage1(:, 1)';
outputStage2MidTemp = fconv(outputStage1MidTemp, midFirHann);
outputStage1SideTemp = outputStage1(:, 2)';
outputStage2SideTemp = fconv(outputStage1SideTemp, sideFirHann);
disp(['Done convolution in ', num2str(toc), ' sec.'])
disp('Cutting convolution edges and converting MS to LR...')
outputStage2Mid = outputStage2MidTemp( (double(fftSize) / 2 + 1) : (size(outputStage2MidTemp, 2) - double(fftSize) / 2 + 1) )';
outputStage2Side = outputStage2SideTemp( (double(fftSize) / 2 + 1) : (size(outputStage2SideTemp, 2) - double(fftSize) / 2 + 1) )';
outputStage2(:, 1) = outputStage2Mid + outputStage2Side;
outputStage2(:, 2) = outputStage2Mid - outputStage2Side;
disp('Partly cleaning memory...')
clearvars -except currentSessionStatus outputStage2 finalAmplitudeCoefficient pathWorkspace fileTarget fileReference pathTemp pathBin fileBrickwallLimiter maxLength timeArea fftSize linLogOversampling peakCompensationSteps targetSamplingFrequency targetDivisions targetPartSize referenceMatchingRms;
disp('Stage 2 completed')
disp('--------------------------------')
% Stage 3 - RMS Correction
disp('Stage 3 - RMS Correction')
currentSessionStatus = 6;
% Check to use limiter or not
usingLimiter = true;
if strcmp(pathBin, '') == true
usingLimiter = false;
disp('Brickwall limiter will not be used')
end
disp('Calculating temp mid channel of matched target...')
outputStage2Mid(:) = (outputStage2(:, 1) + outputStage2(:, 2)) ./ 2;
% Will process outputStage3 in next 'for'
outputStage3 = outputStage2;
for p = 1 : 1 : peakCompensationSteps
% Limiter emulation :)
if p < peakCompensationSteps
outputStage2MidLimited = outputStage2Mid;
disp(['Applying peak compensation #', num2str(p), ': simple algorithm...'])
outputStage2MidLimited( outputStage2MidLimited > 1 ) = 1;
outputStage2MidLimited( outputStage2MidLimited < -1 ) = -1;
end
% Loading brickwall limiter
if p == peakCompensationSteps
if usingLimiter == false
break
end
disp(['Applying peak compensation #', num2str(p), ': complex algorithm...'])
audiowrite([pathTemp, 'tempnolim.wav'], outputStage2Mid, targetSamplingFrequency, 'BitsPerSample', 32);
command = ['"', pathBin, 'mrswatson64.exe" -p "', pathBin, fileBrickwallLimiter, '" -b 2048 -i "', pathTemp, 'tempnolim.wav" -o "', pathTemp, 'templim.wav"'];
[~, ~] = system(command);
[outputStage2MidLimited, ~] = audioread([pathTemp, 'templim.wav']);
delete([pathTemp, 'tempnolim.wav']);
delete([pathTemp, 'templim.wav']);
end
% Some code from Stage 1
outputStage2Rmses = zeros(1, targetDivisions);
for i = 1 : 1 : targetDivisions
outputStage2Rmses(i) = rms(outputStage2MidLimited( (targetPartSize * (i - 1) + 1) : (targetPartSize * i) ));
end
outputStage2AverageRms = rms(outputStage2Rmses);
% May add array preallocation here in future
j = 1;
for i = 1 : 1 : targetDivisions
if outputStage2Rmses(i) >= (outputStage2AverageRms)
outputStage2AnalyzedPartNumbers(j) = i;
j = j + 1;
end
end
outputStage2AnalyzedRmses = zeros(1, size(outputStage2AnalyzedPartNumbers, 2));
for i = 1 : 1 : size(outputStage2AnalyzedPartNumbers, 2)
outputStage2AnalyzedRmses(i) = outputStage2Rmses(outputStage2AnalyzedPartNumbers(i));
end
disp(['Current average RMS value in loudest parts is ', num2str(rms(outputStage2AnalyzedRmses))])
% ---
outputStage2AnalyzedRmses(outputStage2AnalyzedRmses < 0.000001) = 0.000001;
rmsCoefficient = referenceMatchingRms / rms(outputStage2AnalyzedRmses);
if p < peakCompensationSteps
disp(['The RMS correction coefficient #', num2str(p), ' is: ', num2str(rmsCoefficient), '. Modifying amplitudes...'])
outputStage2Mid = outputStage2Mid .* rmsCoefficient;
outputStage3 = outputStage3 .* rmsCoefficient;
end
clearvars outputStage2AnalyzedPartNumbers outputStage2AnalyzedRmses outputStage2MidLimited;
end
disp(['The final RMS correction coefficient is: ', num2str(rmsCoefficient), '. Modifying amplitudes of the target...'])
outputStage3 = outputStage3 .* rmsCoefficient;
disp('Partly cleaning memory...')
clearvars -except currentSessionStatus outputStage3 finalAmplitudeCoefficient usingLimiter pathWorkspace fileTarget fileReference pathTemp pathBin fileBrickwallLimiter maxLength timeArea fftSize linLogOversampling peakCompensationSteps targetSamplingFrequency;
disp('Stage 3 completed')
disp('--------------------------------')
% Final Stage - Export
disp('Final Stage - Export')
currentSessionStatus = 7;
if usingLimiter == true
% If using limiter, let's make 24-bit and 16-bit WAVs
% We need to compensate first 1041 empty samples from brickwall limiter (this magic number is tested only on Voxengo Elephant)
outputSize = size(outputStage3, 1);
% delayCompensation is a number how much samples with 0 we need to add to the out to be like 2048*n+1
delayCompensation = 2048 * (fix(outputSize / 2048) + 1) - outputSize + 1;
disp(['Adding ', num2str(delayCompensation), ' zeroes in the end of modified target to make delay compensation...'])
outputStage3(outputSize + 1 : outputSize + delayCompensation, :) = 0;
% Loading brickwall limiter for final limiting
disp('Limiting modified target...')
audiowrite([pathTemp, 'tempfinnolim.wav'], outputStage3, targetSamplingFrequency, 'BitsPerSample', 32);
command = ['"', pathBin, 'mrswatson64.exe" -p "', pathBin, fileBrickwallLimiter, '" -b 2048 -i "', pathTemp, 'tempfinnolim.wav" -o "', pathTemp, 'tempfinlim.wav"'];
[~, ~] = system(command);
[outputStage4, ~] = audioread([pathTemp, 'tempfinlim.wav']);
delete([pathTemp,'tempfinnolim.wav']);
delete([pathTemp,'tempfinlim.wav']);
% Trimming
disp('Processing delay compensation...')
outputStage4 = outputStage4( (1041 + 1) : (1041 + size(outputStage4, 1) - 2048 - delayCompensation + 1), :);
%clearvars outputStage4;
if finalAmplitudeCoefficient ~= 1
% Modifying with final coefficient
disp(['Modifying amplitudes with final aplitude coefficient: ', num2str(finalAmplitudeCoefficient), '...'])
outputStage4 = outputStage4 .* finalAmplitudeCoefficient;
end
currentSessionStatus = 8;
% Final saving to WAV 24-bit
fileOutput = ['[Matchered 24-bit] ', fileTarget];
disp(['Saving modified target audio in WAV 24bit to ', pathWorkspace, fileOutput, '...'])
audiowrite([pathWorkspace, fileOutput], outputStage4, targetSamplingFrequency, 'BitsPerSample', 24);
disp([fileOutput, ' saved'])
% Final saving to WAV 16-bit
fileOutput = ['[Matchered 16-bit] ', fileTarget];
disp(['Saving modified target audio in WAV 16bit to ', pathWorkspace, fileOutput, '...'])
audiowrite([pathWorkspace, fileOutput], outputStage4, targetSamplingFrequency, 'BitsPerSample', 16);
disp([fileOutput, ' saved'])
currentSessionStatus = 9;
else
% If not using limiter, let's make just 32-bit WAV
if finalAmplitudeCoefficient ~= 1
% Modifying with final coefficient
disp(['Modifying amplitudes with final aplitude coefficient: ', num2str(finalAmplitudeCoefficient), '...'])
outputStage4 = outputStage3 .* finalAmplitudeCoefficient;
else
outputStage4 = outputStage3;
end
currentSessionStatus = 8;
% Final saving to WAV 32-bit
fileOutput = ['[Matchered 32-bit NO LIMITING] ', fileTarget];
disp(['Saving modified target audio in WAV 32bit to ', pathWorkspace, fileOutput, '...'])
audiowrite([pathWorkspace, fileOutput], outputStage4, targetSamplingFrequency, 'BitsPerSample', 32);
disp([fileOutput, ' saved'])
currentSessionStatus = 9;
end
currentSessionStatus = 9;
disp('Final stage completed')
end
end
end
end
% Cleaning for safe...
clearvars -except currentSessionStatus
disp('----------------------------------------------------------------')
disp('Stopping DSP Module...')
end