From a7ebf4cdb176f52f5004986b9d8da6fa97b1241f Mon Sep 17 00:00:00 2001 From: Sven Franz Date: Fri, 29 May 2020 13:58:21 +0200 Subject: [PATCH] update to Matlab 2018b + minor changes (try/catch) --- +NoiseSynthesis/+external/DeCrackleNoise.m | 6 +++++- +NoiseSynthesis/+external/STFT.m | 4 ++-- .../+external/computeBandMixingMatrix.m | 1 + .../NoiseAnalysisSynthesis.m | 6 +++++- +NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m | 1 + .../private/analyzeAmplitudeDistribution.m | 13 +++++++++++-- .../private/applyAmplitudeDistribution.m | 1 + .../private/computeLevelFluctuations.m | 2 ++ .../private/computeMixingMatrix.m | 6 +++++- .../private/estimateClickBandwidth.m | 4 ++-- .../private/manualCoherence.m | 15 +++++++++++++++ .../private/shapeAmplitudeDistr.m | 11 +++++++++-- +NoiseSynthesis/ErrorMeasures.m | 4 ++-- +NoiseSynthesis/ModelParametersSet.m | 4 +++- 14 files changed, 64 insertions(+), 14 deletions(-) create mode 100644 +NoiseSynthesis/@NoiseAnalysisSynthesis/private/manualCoherence.m diff --git a/+NoiseSynthesis/+external/DeCrackleNoise.m b/+NoiseSynthesis/+external/DeCrackleNoise.m index 3613a96..c3368ce 100644 --- a/+NoiseSynthesis/+external/DeCrackleNoise.m +++ b/+NoiseSynthesis/+external/DeCrackleNoise.m @@ -26,7 +26,11 @@ end stAlgo = DeCrackle(fs,ThreshDeCrackle); stAlgo.init(); -[vDeCrackled] = stAlgo.process(vSignal); +try + [vDeCrackled] = stAlgo.process(vSignal); +catch + vDeCrackled = vSignal; +end diff --git a/+NoiseSynthesis/+external/STFT.m b/+NoiseSynthesis/+external/STFT.m index f8818b7..7eeb75d 100644 --- a/+NoiseSynthesis/+external/STFT.m +++ b/+NoiseSynthesis/+external/STFT.m @@ -106,14 +106,14 @@ if nargout, varargout = {mSpectrogram, vFreq, vTime}; else - mPSD = (mSpectrogram .* conj(mSpectrogram)) / (vWindow'*vWindow) / fs; + mPSD = (mSpectrogram .* (mSpectrogram)) / (vWindow'*vWindow) / fs; if strcmpi(szSpecType,'single'), mPSD = mPSD(1:iFShalf,:); mPSD(2:end-1) = 2*mPSD(2:end-1); end - surf(vTime,vFreq,10*log10(mPSD + eps^2),'edgecolor','none','AmbientStrength',0.5); + surf(vTime,vFreq,10*log10(abs(mPSD) + eps^2),'edgecolor','none','AmbientStrength',0.5); lightangle(210, 45); axis([vTime(1),vTime(end),vFreq(1),vFreq(end)]); diff --git a/+NoiseSynthesis/+external/computeBandMixingMatrix.m b/+NoiseSynthesis/+external/computeBandMixingMatrix.m index 34d3a10..834cc40 100644 --- a/+NoiseSynthesis/+external/computeBandMixingMatrix.m +++ b/+NoiseSynthesis/+external/computeBandMixingMatrix.m @@ -18,6 +18,7 @@ function [mMixingMatrix] = computeBandMixingMatrix(mGamma) % +mGamma(isnan(mGamma)) = eps; [mEigVecs,mEigVals] = eig(mGamma); mMixingMatrix = (sign(mEigVals) .* sqrt(abs(mEigVals))) * mEigVecs'; diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/NoiseAnalysisSynthesis.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/NoiseAnalysisSynthesis.m index 41b9287..c728660 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/NoiseAnalysisSynthesis.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/NoiseAnalysisSynthesis.m @@ -37,7 +37,7 @@ properties ( Access = private ) vOriginalAnalysisSignal; % Raw analysis signal vBeforeDeCrackling; - mBands; % Frequency bands +% mBands; % Frequency bands mLevelCurvesDecorr; % Decorrelated level fluctuations for all bands blocklen = 1024; % Block length for STFT @@ -95,6 +95,8 @@ properties ( Access = public ) bEstimateClickSpec = true; % Bool whether to estimate the click cutoff frequency in the analysis. Dependent on bApplyClicks [default: true] bDeClick = true; % Bool whether to declick the analysis file [default: true] + + mBands; % Frequency bands end properties ( Access = private, Dependent ) @@ -235,6 +237,8 @@ methods hCohereFun = @(Freq,dist,theta,vPSD) binaural2d(self,dist,Freq); case 'binaural3d', hCohereFun = @(Freq,dist,theta,vPSD) binaural3d(self,dist,Freq); + case 'manual', + hCohereFun = @(Freq,dist,theta,vPSD) manualCoherence(self,Freq,dist); otherwise, warning(sprintf('Coherence model not recognized. Switched to default (''%s'')...',... self.ModelParameters.CohereModel)); %#ok diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m index 9b7e8de..3c2364a 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m @@ -36,6 +36,7 @@ if self.bHPFilterAnalysis, self.AnalysisSignal = filter(b,a,self.AnalysisSignal); end + % estimate the amplitude distribution showMsg(self,'Analyzing Amplitude Distribution'); analyzeAmplitudeDistribution(self); diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/analyzeAmplitudeDistribution.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/analyzeAmplitudeDistribution.m index f69af56..be528e5 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/analyzeAmplitudeDistribution.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/analyzeAmplitudeDistribution.m @@ -122,10 +122,15 @@ switch self.ModelParameters.AmplitudeModel, if isempty(perclower), perclower = 0.05; else + a = perclower; perclower = perclower(... perclower >= percRangeLower(1) & ... perclower <= percRangeLower(2)); - perclower = perclower(1); + if ~isempty(perclower), + perclower = perclower(1); + else + perclower = 0.05; + end end if isempty(percupper), percupper = 0.95; @@ -133,7 +138,11 @@ switch self.ModelParameters.AmplitudeModel, percupper = percupper(... percupper >= percRangeUpper(1) & ... percupper <= percRangeUpper(2)); - percupper = percupper(end); + if ~isempty(percupper) + percupper = percupper(end); + else + percupper = 0.95; + end end % warning('off','stats:paretotails:ConvergedToBoundaryLower'); diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/applyAmplitudeDistribution.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/applyAmplitudeDistribution.m index a908c6b..d580423 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/applyAmplitudeDistribution.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/applyAmplitudeDistribution.m @@ -14,6 +14,7 @@ import NoiseSynthesis.external.* if self.bApplyAmplitudeDistr, showMsg(self,'Applying desired Amplitude Distribution'); + self.SensorSignals(isnan(self.SensorSignals)) = 0; for iSignal = 1:self.NumSensorSignals, self.SensorSignals(:, iSignal) = ... shapeAmplitudeDistr(self, self.SensorSignals(:, iSignal)); diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeLevelFluctuations.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeLevelFluctuations.m index 48104d0..757bc9a 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeLevelFluctuations.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeLevelFluctuations.m @@ -22,6 +22,8 @@ vIdxNormalize = ... round(0.05 * self.numBlocks) : ... round(0.95 * self.numBlocks); +vIdxNormalize(vIdxNormalize <= 0) = []; + self.mLevelCurves = zeros(self.lenLevelCurve,self.numBands); for aaBand = 1:self.numBands, vCurrBandSignal = self.mBands(aaBand,:).'; diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeMixingMatrix.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeMixingMatrix.m index aaaa486..c2d2f25 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeMixingMatrix.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/computeMixingMatrix.m @@ -32,7 +32,11 @@ for ppSensor = 1:self.NumSensorSignals, end end -[mEigVecs,mEigVals] = eig(mGamma); +try + [mEigVecs,mEigVals] = eig(mGamma); +catch + disp("TEST"); +end mEigVals(abs(mEigVals) <= eps) = 0; % prevent numerical issues mMixingMatrix = sqrt(mEigVals) * mEigVecs'; diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/estimateClickBandwidth.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/estimateClickBandwidth.m index 718883e..9857256 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/estimateClickBandwidth.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/estimateClickBandwidth.m @@ -11,8 +11,8 @@ function [] = estimateClickBandwidth(self,vClicks) blocklenSec = 32e-3; overlapRatio = 0.5; -params = STFTparams(blocklenSec,overlapRatio,self.Fs); -mSTFT = STFT(vClicks,params); +params = NoiseSynthesis.STFTparams(blocklenSec,overlapRatio,self.Fs); +mSTFT = NoiseSynthesis.external.STFT(vClicks,params); vFreq = linspace(0,self.Fs/2,params.NFFT/2+1); diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/manualCoherence.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/manualCoherence.m new file mode 100644 index 0000000..7903bf2 --- /dev/null +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/manualCoherence.m @@ -0,0 +1,15 @@ +function [vCohere] = manualCoherence(self,vFreq, dist) + +if ~isempty(self.ModelParameters.manualCohereModel) + vFreq = vFreq(:); + blocklen = (size(self.ModelParameters.MeanPSD, 1) - 1) * 2; + + vFreqIdx = floor(vFreq * blocklen / self.Fs + 1); + vCohere = self.ModelParameters.manualCohereModel(vFreqIdx); + if dist == 0 || isnan(vCohere) + vCohere = ones(size(vCohere)); + end +else + vCohere = 1; +end + diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/shapeAmplitudeDistr.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/shapeAmplitudeDistr.m index 81e8350..8414836 100644 --- a/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/shapeAmplitudeDistr.m +++ b/+NoiseSynthesis/@NoiseAnalysisSynthesis/private/shapeAmplitudeDistr.m @@ -14,12 +14,20 @@ import NoiseSynthesis.external.* % 'inverse method' or 'percentile transformation method' stated % in Papoulis. +% if ~isempty(isnan(vSigIn)) +% vSigOut = vSigIn .* 0; +% return +% end [vCDFsignal,vXSignal] = ecdf(vSigIn); vXSignal = vXSignal(2:end); vCDFsignal = (vCDFsignal(1:end-1) + vCDFsignal(2:end))/2; -vUniformlyDistrNoise = interp1(vXSignal,vCDFsignal,vSigIn,'linear','extrap'); +try + vUniformlyDistrNoise = interp1(vXSignal,vCDFsignal,vSigIn,'linear','extrap'); +catch + disp("TEST"); +end numPoints = 1000; switch self.ModelParameters.AmplitudeModel, @@ -77,5 +85,4 @@ vSigOut = interp1(vCDF,vQuantiles,... 'linear','extrap'); - % End of file: shapeAmplitudeDistr.m diff --git a/+NoiseSynthesis/ErrorMeasures.m b/+NoiseSynthesis/ErrorMeasures.m index 5f98084..21e204e 100644 --- a/+NoiseSynthesis/ErrorMeasures.m +++ b/+NoiseSynthesis/ErrorMeasures.m @@ -178,7 +178,7 @@ classdef ErrorMeasures < handle case 'cosh', % from % Gray and Markel, Distance Measures for Speech Processing - error = distchpf(vPSDReference(:).',vPSDSynthesis(:).'); + error = NoiseSynthesis.external.distchpf(vPSDReference(:).',vPSDSynthesis(:).'); case 'logMSE', error = mean( ... (log10(abs(vPSDReference)) - log10(abs(vPSDSynthesis))).^2 ... @@ -241,7 +241,7 @@ classdef ErrorMeasures < handle error = KullbackLeibler(quantilesP, pdfP, pdfQ); case 'Bhattacharyya', - error = bhattacharyya(pdfP(:), pdfQ(:)); + error = NoiseSynthesis.external.bhattacharyya(pdfP(:), pdfQ(:)); case 'Hellinger', error = sqrt(1 - bhattacharyya(pdfP, pdfQ)); diff --git a/+NoiseSynthesis/ModelParametersSet.m b/+NoiseSynthesis/ModelParametersSet.m index ebaf091..872921e 100644 --- a/+NoiseSynthesis/ModelParametersSet.m +++ b/+NoiseSynthesis/ModelParametersSet.m @@ -43,7 +43,8 @@ properties ( Transient, Hidden ) 'spherical',... 'anisotropic',... 'binaural2d',... - 'binaural3d'}); + 'binaural3d',... + 'manual'}); end properties ( Access = public ) @@ -71,6 +72,7 @@ properties ( Access = public ) CDF; % Dependent variable of ECDF or holds information for other models CohereModel = 'cylindrical'; % Spatial coherence model + manualCohereModel = []; % Sensor position(s) in meters SensorPositions = [ 0 0.17 -- 2.18.1