diff --git a/+NoiseSynthesis/+external/DeCrackleNoise.m b/+NoiseSynthesis/+external/DeCrackleNoise.m
index 3613a9694a1293b71f6f60bb7104042974bcd9f4..c3368cecc827b1e2ced47233a35421bb5900ff86 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 f8818b77ca8d629df13fcdd622181eb767ba0063..7eeb75d6223bb029e6ecbbf2a2d1da7e610c7d10 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 34d3a1086410bf3069582811129a3c78b98e758f..834cc406d7a1349f98aee432280fb2ffdddfdd2b 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 41b9287d41e66b77a8321e78fc12802088d571df..c728660be4ae0fcf9b64688a49f6781dc9ad0f68 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<SPWRN>
diff --git a/+NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m b/+NoiseSynthesis/@NoiseAnalysisSynthesis/analyze.m
index 9b7e8de236e9ae56666964317dd8119f2e97cec0..3c2364abec5a205a75e33a74ec3ddc5db347631c 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 f69af56374d8f4d150df1317f9eea625386aed12..be528e589bec4ea0bd2806537def8c79efb41ce1 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 a908c6bd80e691e70b302628cb12cea63fda052e..d5804230b5b1e3b5db645f48405dd06528723aa1 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 48104d0cb06bebc9403add5f5a717aa599fcf34e..757bc9aadd374b707cbcd912b529a013f38b722a 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 aaaa4865d516e3ca3e06313cfc0b3fe6bd9ed16e..c2d2f25489a95557ace79876dc703bd78b399e49 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 718883eb2580adc6d87044c00eca0e9f7be397f6..9857256d4aed78d8c5e67a5d361f80512facba41 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 0000000000000000000000000000000000000000..7903bf298d227c63537bb6b614639aa9dd085d0b
--- /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 81e83508d9d174aec9ff1af00749e90110299d1b..8414836f610507260951166e95c17e79ec32ea3a 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 5f980840b93894ed3e820ebfceaf8df945963415..21e204efdbac38ee806c1a5f9ded213a596ed758 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 ebaf091fbaa7f04d78fc95b9b96aa748359d9dae..872921eca63a60e1ff0140e735da289de73a2c0f 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