diff --git a/applications/utilities/postProcessing/noise/noise.C b/applications/utilities/postProcessing/noise/noise.C
index 5c3ad2a139..cf5794075f 100644
--- a/applications/utilities/postProcessing/noise/noise.C
+++ b/applications/utilities/postProcessing/noise/noise.C
@@ -107,7 +107,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
 
-    word dictName("noiseDict");
+    fileName dictName(runTime.system()/"noiseDict");
     if (args.optionFound("dict"))
     {
         dictName = args["dict"];
@@ -118,7 +118,6 @@ int main(int argc, char *argv[])
         IOobject
         (
             dictName,
-            runTime.system(),
             runTime,
             IOobject::MUST_READ
         )
diff --git a/applications/utilities/postProcessing/noise/noiseDict b/applications/utilities/postProcessing/noise/noiseDict
new file mode 100644
index 0000000000..e983cabc12
--- /dev/null
+++ b/applications/utilities/postProcessing/noise/noiseDict
@@ -0,0 +1,136 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  plus                                  |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      noiseDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+noiseModel      surfaceNoise;
+
+surfaceNoiseCoeffs
+{
+    windowModel     Hanning;
+
+    HanningCoeffs
+    {
+        // Window overlap percentage
+        overlapPercent  50;
+        symmetric       yes;
+        extended        yes;
+
+        // Optional number of windows, default = all available
+        // nWindow         1;
+    }
+
+/*
+    windowModel     uniform;
+
+    uniformCoeffs
+    {
+        // Window overlap percentage
+        overlapPercent  50;
+
+        value           1;
+
+        // Optional number of windows, default = all available
+        // nWindow         1;
+    }
+*/
+
+
+    // Input file
+    inputFile   "postProcessing/faceSource1/surface/patch_motorBike_rider-helmet%65/patch_motorBike_rider-helmet%65.case";
+
+    // Surface reader
+    reader      ensight;
+
+    // Surface writer
+    writer      ensight;
+
+    // Collate times for ensight output - ensures geometry is only written once
+    writeOptions
+    {
+        ensight
+        {
+            collateTimes 1;
+        }
+    }
+
+    // Reference density (to convert from kinematic to static pressure)
+    rhoRef          1.205;
+
+    // Number of samples in sampling window
+    // Must be a power of 2, default = 2^16 (=65536)
+    N               4096; // 8192; // 4096;
+
+    // Lower frequency limit, default = 25Hz
+    //fl              25;
+
+    // Upper frequency limit, default = 10kHz
+    fu              15000;
+
+    // Start time, default = 0s
+    //startTime       0;
+
+    // Write interval for FFT data, default = 1
+//    fftWriteInterval 100;
+}
+
+pointNoiseCoeffs
+{
+    csvFileData
+    {
+        fileName        "pressureData";
+        nHeaderLine     1;
+        refColumn       0;
+        componentColumns (1);
+        separator       " ";
+        mergeSeparators yes;
+    }
+
+    HanningCoeffs
+    {
+        // Window overlap percentage
+        overlapPercent  50;
+        symmetric       yes;
+        extended        yes;
+
+        // Optional number of windows, default = all available
+        //nWindow         5;
+    }
+
+    // Graph format, default = raw
+    graphFormat     raw;
+
+    // Reference density (to convert from kinematic to static pressure)
+    rhoRef          1.2;
+
+    // Number of samples in sampling window
+    // Must be a power of 2, default = 2^16 (=65536)
+    N               4096;
+
+    // Lower frequency limit, default = 25Hz
+    //fl              25;
+
+    // Upper frequency limit, default = 10kHz
+    //fu              10000;
+
+    // Start time, default = 0s
+    //startTime       0;
+
+    // Write interval for FFT data, default = 1
+    fftWriteInterval 100;
+}
+
+
+// ************************************************************************* //
diff --git a/src/randomProcesses/Make/files b/src/randomProcesses/Make/files
index a16bde251c..11138bd932 100644
--- a/src/randomProcesses/Make/files
+++ b/src/randomProcesses/Make/files
@@ -3,7 +3,6 @@ $(Kmesh)/Kmesh.C
 
 fft = fft
 $(fft)/fft.C
-$(fft)/fftRenumber.C
 $(fft)/calcEk.C
 $(fft)/kShellIntegration.C
 
@@ -26,6 +25,7 @@ windowModels = windowModels
 $(windowModels)/windowModel/windowModel.C
 $(windowModels)/windowModel/windowModelNew.C
 $(windowModels)/Hanning/Hanning.C
+$(windowModels)/uniform/uniform.C
 
 
 LIB = $(FOAM_LIBBIN)/librandomProcesses
diff --git a/src/randomProcesses/Make/options b/src/randomProcesses/Make/options
index 86d0d7baee..e1eddd51d2 100644
--- a/src/randomProcesses/Make/options
+++ b/src/randomProcesses/Make/options
@@ -1,9 +1,11 @@
 EXE_INC = \
+    -I$(FFTW_ARCH_PATH)/include \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
+    -L$(FFTW_ARCH_PATH)/lib$(WM_COMPILER_LIB_ARCH) -lfftw3 \
     -lfiniteVolume \
     -lsampling \
     -lsurfMesh
diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.C b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
index 0d36c47d16..abaf199514 100644
--- a/src/randomProcesses/noise/noiseFFT/noiseFFT.C
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.C
@@ -26,9 +26,10 @@ License
 #include "noiseFFT.H"
 #include "IFstream.H"
 #include "DynamicList.H"
-#include "fft.H"
 #include "SubField.H"
 #include "mathematicalConstants.H"
+#include "HashSet.H"
+#include "fft.H"
 
 using namespace Foam::constant;
 
@@ -44,7 +45,7 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
 )
 {
     tmp<scalarField> tf(new scalarField(N/2, 0));
-    scalarField& f = tf();
+    scalarField& f = tf.ref();
 
     scalar deltaf = 1.0/(N*deltaT);
     forAll(f, i)
@@ -56,70 +57,80 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::frequencies
 }
 
 
-void Foam::noiseFFT::octaveFrequenciesIDs
+Foam::tmp<Foam::scalarField> Foam::noiseFFT::PSD(const scalarField& PSDf)
+{
+    return 10*log10(PSDf/sqr(p0));
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::noiseFFT::SPL(const scalarField& Prms2)
+{
+    return 10*log10(Prms2/sqr(p0));
+}
+
+
+void Foam::noiseFFT::octaveBandInfo
 (
     const scalarField& f,
     const scalar fLower,
     const scalar fUpper,
     const scalar octave,
-    labelList& freqBandIDs
+    labelList& fBandIDs,
+    scalarField& fCentre
 )
 {
     // Set the indices of to the lower frequency bands for the input frequency
     // range. Ensure that the centre frequency passes though 1000 Hz
 
     // Low frequency bound given by:
-    //     fLow = f0*(2^(bandI/octave/2))
-    // Centre frequency given by:
-    //     fCentre = f0*(2^(bandI/octave))
+    //     fLow = f0*(2^(0.5*bandI/octave))
 
-    scalar f0 = 1000;
-    scalar minFrequency = max(fLower, min(f));
+    // Initial (lowest centre frequency)
+    scalar fTest = 15.625;
 
-    // Lower frequency band limit
-    label band0Low = ceil(2*octave*log(minFrequency/f0)/log(2.0));
+    const scalar fRatio = pow(2, 1.0/octave);
+    const scalar fRatioL2C = pow(2, 0.5/octave);
 
-    // Centre frequency band limit
-    //label band0Centre = ceil(octave*log(fLower/f0)/log(2.0));
+    // IDs of band IDs
+    labelHashSet bandIDs(f.size());
 
-    scalar fLowerBand = f0*pow(2, band0Low/octave/2);
-    scalar fRatio = pow(2, 1.0/octave);
+    // Centre frequencies
+    DynamicList<scalar> fc;
+
+    // Convert to lower band limit
+    fTest /= fRatioL2C;
 
-    bool complete = false;
-    DynamicList<label> bandIDs(f.size());
     forAll(f, i)
     {
-        while (f[i] >= fLowerBand)
+        if (f[i] >= fTest)
         {
-            bandIDs.append(i);
-            fLowerBand *= fRatio;
-
-            if (fLowerBand > fUpper)
+            // Advance band if appropriate
+            while (f[i] > fTest)
+            {
+                fTest *= fRatio;
+            }
+            fTest /= fRatio;
+
+            if (bandIDs.insert(i))
+            {
+                // Also store (next) centre frequency
+                fc.append(fTest*fRatioL2C);
+            }
+            fTest *= fRatio;
+
+            if (fTest > fUpper)
             {
-                complete = true;
                 break;
             }
         }
-
-        if (complete) break;
     }
 
-    freqBandIDs.transfer(bandIDs);
-}
+    fBandIDs = bandIDs.toc();
 
+    // Remove the last centre frequency (beyond upper frequency limit)
+    fc.remove();
 
-Foam::tmp<Foam::scalarField> Foam::noiseFFT::octaveFrequencies
-(
-    const scalarField& f,
-    const scalar fLower,
-    const scalar fUpper,
-    const scalar octave
-)
-{
-    labelList freqBandIDs;
-    octaveFrequenciesIDs(f, fLower, fUpper, octave, freqBandIDs);
-    tmp<scalarField> tf(new scalarField(f, freqBandIDs));
-    return tf;
+    fCentre.transfer(fc);
 }
 
 
@@ -133,7 +144,10 @@ Foam::noiseFFT::noiseFFT
 :
     scalarField(pressure),
     deltaT_(deltaT)
-{}
+{
+    scalarField& p = *this;
+    p -= average(p);
+}
 
 
 Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
@@ -191,6 +205,9 @@ Foam::noiseFFT::noiseFFT(const fileName& pFileName, const label skip)
     deltaT_ = (T1 - T0)/pData.size();
 
     this->transfer(pData);
+
+    scalarField& p = *this;
+    p -= average(p);
 }
 
 
@@ -220,6 +237,8 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
     const tmp<scalarField>& tpn
 ) const
 {
+    // Calculate the 2-sided fft
+    // Note: result not scaled
     tmp<scalarField> tPn2
     (
         mag
@@ -234,17 +253,15 @@ Foam::tmp<Foam::scalarField> Foam::noiseFFT::Pf
 
     tpn.clear();
 
+    // Only storing the positive half
+    // Note: storing (N/2) values, DC component at position (0)
     tmp<scalarField> tPn
     (
         new scalarField
         (
-            scalarField::subField(tPn2(), tPn2().size()/2)
+            scalarField::subField(tPn2(), tPn2().size()/2 + 1)
         )
     );
-    scalarField& Pn = tPn();
-
-    Pn *= 2.0/sqrt(scalar(tPn2().size()));
-    Pn[0] /= 2.0;
 
     return tPn;
 }
@@ -255,7 +272,7 @@ Foam::graph Foam::noiseFFT::meanPf(const windowModel& window) const
     const label N = window.nSamples();
     const label nWindow = window.nWindow();
 
-    scalarField meanPf(N/2, 0.0);
+    scalarField meanPf(N/2 + 1, 0.0);
 
     for (label windowI = 0; windowI < nWindow; ++windowI)
     {
@@ -287,14 +304,13 @@ Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
     const label N = window.nSamples();
     const label nWindow = window.nWindow();
 
-    scalarField RMSMeanPf(N/2, 0.0);
-
+    scalarField RMSMeanPf(N/2 + 1, 0.0);
     for (label windowI = 0; windowI < nWindow; ++windowI)
     {
         RMSMeanPf += sqr(Pf(window.apply<scalar>(*this, windowI)));
     }
 
-    RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow));
+    RMSMeanPf = sqrt(RMSMeanPf/scalar(nWindow))/scalar(N);
 
     scalar deltaf = 1.0/(N*deltaT_);
     scalarField f(RMSMeanPf.size());
@@ -305,9 +321,9 @@ Foam::graph Foam::noiseFFT::RMSmeanPf(const windowModel& window) const
 
     return graph
     (
-        "P(f)",
+        "Prms(f)",
         "f [Hz]",
-        "P(f) [Pa]",
+        "Prms(f) [Pa]",
         f,
         RMSMeanPf
     );
@@ -319,19 +335,30 @@ Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const
     const label N = window.nSamples();
     const label nWindow = window.nWindow();
 
-    scalarField psd(N/2, 0.0);
+    scalarField psd(N/2 + 1, 0.0);
 
     for (label windowI = 0; windowI < nWindow; ++windowI)
     {
-        psd += 0.5*sqr(Pf(window.apply<scalar>(*this, windowI)));
+        psd += sqr(Pf(window.apply<scalar>(*this, windowI)));
     }
 
     scalar deltaf = 1.0/(N*deltaT_);
+    scalar fs =  1.0/deltaT_;
+    psd /= scalar(nWindow)*fs*N;
 
-    psd /= nWindow*deltaf;
+    // Scaling due to use of 1-sided FFT
+    // Note: do not scale DC component
+    psd *= 2;
+    psd.first() /= 2;
+    psd.last() /= 2;
 
     scalarField f(psd.size());
 
+    if (0) // if (debug)
+    {
+        Pout<< "Average PSD: " << average(psd) << endl;
+    }
+
     forAll(f, i)
     {
         f[i] = i*deltaf;
@@ -348,159 +375,87 @@ Foam::graph Foam::noiseFFT::PSDf(const windowModel& window) const
 }
 
 
-Foam::graph Foam::noiseFFT::PSD(const graph& gPSD) const
+Foam::graph Foam::noiseFFT::PSD(const graph& gPSDf) const
 {
     return graph
     (
-        "PSD(dB)",
+        "PSD(f)",
         "f [Hz]",
-        "PSD_dB(f) [dB]",
-        gPSD.x(),
-        10*log10(gPSD.y()/sqr(p0))
+        "PSD_dB(f) [dB_Hz]",
+        gPSDf.x(),
+        10*log10(gPSDf.y()/sqr(p0))
     );
 }
 
 
-Foam::graph Foam::noiseFFT::Lf(const graph& gPf) const
-{
-    return graph
-    (
-        "L(f)",
-        "f [Hz]",
-        "L(f) [dB]",
-        gPf.x(),
-        20*log10(gPf.y()/p0)
-    );
-}
-
-
-Foam::graph Foam::noiseFFT::Ldelta
+Foam::graph Foam::noiseFFT::octaves
 (
-    const graph& gLf,
-    const labelList& freqBandIDs
+    const graph& g,
+    const labelList& freqBandIDs,
+    bool integrate
 ) const
 {
     if (freqBandIDs.size() < 2)
     {
         WarningInFunction
             << "Octave frequency bands are not defined "
-            << "- skipping Ldelta calculation"
+            << "- skipping octaves calculation"
             << endl;
 
         return graph
         (
-            "Ldelta",
-            "fm [Hz]",
-            "Ldelta [dB]",
+            "octave",
+            "x",
+            "y",
             scalarField(),
             scalarField()
         );
     }
 
-    const scalarField& f = gLf.x();
-    const scalarField& Lf = gLf.y();
+    const scalarField& f = g.x();
+    const scalarField& data = g.y();
 
-    scalarField ldelta(freqBandIDs.size() - 1, 0.0);
+    scalarField octData(freqBandIDs.size() - 1, 0.0);
     scalarField fm(freqBandIDs.size() -1, 0.0);
 
     for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
     {
-        label f0 = freqBandIDs[bandI];
-        label f1 = freqBandIDs[bandI+1];
-        fm[bandI] = f[f0];
+        label fb0 = freqBandIDs[bandI];
+        label fb1 = freqBandIDs[bandI+1];
+        fm[bandI] = f[fb0];
 
-        if (f0 == f1) continue;
+        if (fb0 == fb1) continue;
 
-        for (label freqI = f0; freqI < f1; freqI++)
+        if (integrate)
         {
-            ldelta[bandI] += pow(10, Lf[freqI]/10.0);
+            for (label freqI = fb0; freqI < fb1; freqI++)
+            {
+                label f0 = f[freqI];
+                label f1 = f[freqI + 1];
+                scalar dataAve = 0.5*(data[freqI] + data[freqI + 1]);
+                octData[bandI] += dataAve*(f1 - f0);
+            }
+        }
+        else
+        {
+            for (label freqI = fb0; freqI < fb1; freqI++)
+            {
+                octData[bandI] += data[freqI];
+            }
         }
-        ldelta[bandI] = 10*log10(ldelta[bandI]);
     }
 
     return graph
     (
-        "Ldelta",
+        "octaves(f)",
         "fm [Hz]",
-        "Ldelta [dB]",
+        "octave data",
         fm,
-        ldelta
+        octData
     );
 }
 
 
-Foam::graph Foam::noiseFFT::Pdelta
-(
-    const graph& gPf,
-    const labelList& freqBandIDs
-) const
-{
-    if (freqBandIDs.size() < 2)
-    {
-        WarningInFunction
-            << "Octave frequency bands are not defined "
-            << "- skipping Pdelta calculation"
-            << endl;
-
-        return graph
-        (
-            "Pdelta",
-            "fm [Hz]",
-            "Pdelta [dB]",
-            scalarField(),
-            scalarField()
-        );
-    }
-
-    const scalarField& f = gPf.x();
-    const scalarField& Pf = gPf.y();
-
-    scalarField pdelta(freqBandIDs.size() - 1, 0.0);
-    scalarField fm(pdelta.size());
-
-    for (label bandI = 0; bandI < freqBandIDs.size() - 1; bandI++)
-    {
-        label f0 = freqBandIDs[bandI];
-        label f1 = freqBandIDs[bandI+1];
-        fm[bandI] = f[f0];
-
-        if (f0 == f1) continue;
-
-        for (label freqI = f0; freqI < f1; freqI++)
-        {
-            pdelta[bandI] += sqr(Pf[freqI]);
-        }
-        pdelta[bandI] = sqrt((2.0/3.0)*pdelta[bandI]);
-    }
-
-    return graph
-    (
-        "Pdelta",
-        "fm [Hz]",
-        "Pdelta [dB]",
-        fm,
-        pdelta
-    );
-}
-
-
-Foam::scalar Foam::noiseFFT::Lsum(const graph& gLf) const
-{
-    const scalarField& Lf = gLf.y();
-
-    scalar lsum = 0.0;
-
-    forAll(Lf, i)
-    {
-        lsum += pow(10, Lf[i]/10.0);
-    }
-
-    lsum = 10*log10(lsum);
-
-    return lsum;
-}
-
-
 Foam::scalar Foam::noiseFFT::dbToPa(const scalar db) const
 {
     return p0*pow(10.0, db/20.0);
diff --git a/src/randomProcesses/noise/noiseFFT/noiseFFT.H b/src/randomProcesses/noise/noiseFFT/noiseFFT.H
index ab580cbcdf..6439a2f429 100644
--- a/src/randomProcesses/noise/noiseFFT/noiseFFT.H
+++ b/src/randomProcesses/noise/noiseFFT/noiseFFT.H
@@ -32,12 +32,12 @@ Description
     - meanPf: multi-window mean fft
     - RMSmeanPf: multi-window RMS mean fft
     - PSDf: multi-window power spectral density (PSD) in frequency domain
-    - PSD: multi-window power spectral density (PSD) in dB
-    - Lf: narrow-band pressure-fluctuation level (PFL) in frequency domain
+    - PSD: power spectral density in dB/Hz
+    - SPL: sound pressure level in dB
 
     Octave-based data:
-    - Ldelta: PFL spectrum
-    - Pdelta: pressure spectrum
+    - PSD spectrum
+    - SPL spectrum
 
 SourceFiles
     noiseFFT.C
@@ -72,6 +72,17 @@ class noiseFFT
         //- Time spacing of the raw data
         scalar deltaT_;
 
+    struct octaveBandInfo
+    {
+        label octave;
+
+        // IDs of bin boundaries in pressure data
+        labelList binIDs;
+
+        // Centre frequencies for each bin
+        scalarField centreFreq;
+    };
+
 
 public:
 
@@ -104,27 +115,27 @@ public:
             const scalar deltaT
         );
 
+        //- Return the PSD [dB/Hz]
+        //  Input PSD in [Pa^2/Hz]
+        static tmp<scalarField> PSD(const scalarField& PSDf);
+
+        //- Return the SPL [dB]
+        //  Input P(rms)^2 in [Pa^2]
+        static tmp<scalarField> SPL(const scalarField& Prms2);
+
         //- Return a list of the frequency indices wrt f field that
         //  correspond to the bands limits for a given octave
-        static void octaveFrequenciesIDs
+        static void octaveBandInfo
         (
             const scalarField& f,
             const scalar fLower,
             const scalar fUpper,
             const scalar octave,
-            labelList& freqBandIDs
+            labelList& fBandIDs,
+            scalarField& fCentre
         );
 
-        //- Return the 1/octave octave frequency bounds
-        static tmp<scalarField> octaveFrequencies
-        (
-            const scalarField& f,
-            const scalar fLower,
-            const scalar fUpper,
-            const scalar octave
-        );
-
-        //- Return the graph of p(t)
+        //- Return the graph of pressure as a function of time
         graph pt() const;
 
         //- Return the fft of the given pressure data
@@ -141,32 +152,18 @@ public:
         //  pressure data [Pa^2/Hz]
         graph PSDf(const windowModel& window) const;
 
-        //- Return the PSD [dB]
-        graph PSD(const graph& gPSD) const;
+        //- Return the PSD [dB/Hz]
+        //  Takes PSD in [Pa^2/Hz]
+        graph PSD(const graph& gPSDf) const;
 
-        //- Return the narrow-band PFL (pressure-fluctuation level)
-        //  spectrum [dB]
-        graph Lf(const graph& gPf) const;
-
-        //- Return the octave-band PFL spectrum starting at octave
-        //  frequencies given by the supplied frequency bands [dB]
-        graph Ldelta
+        //- Generate octave data
+        graph octaves
         (
-            const graph& gLf,
-            const labelList& freqBandIDs
+            const graph& g,
+            const labelList& freqBandIDs,
+            bool integrate
         ) const;
 
-        //- Return the octave-band pressure spectrum at octave
-        //  frequencies given by the supplied frequency bands [dB]
-        graph Pdelta
-        (
-            const graph& gLf, 
-            const labelList& freqBandIDs
-        ) const;
-
-        //- Return the total PFL as the sum of Lf over all frequencies
-        scalar Lsum(const graph& gLf) const;
-
         //- Convert the db into Pa
         scalar dbToPa(const scalar db) const;
 
diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C
index c32c4ec21d..145f992725 100644
--- a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C
+++ b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2016 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,12 +97,13 @@ Foam::label Foam::noiseModel::findStartTimeIndex
 Foam::noiseModel::noiseModel(const dictionary& dict)
 :
     dict_(dict),
-    pRef_(dict.lookupOrDefault("pRef", 0)),
-    nSamples_(dict.lookupOrDefault("N", 65536)),
-    fLower_(dict.lookupOrDefault("fl", 25)),
-    fUpper_(dict.lookupOrDefault("fu", 10000)),
-    startTime_(dict.lookupOrDefault("startTime", 0)),
-    windowModelPtr_(windowModel::New(dict, nSamples_))
+    rhoRef_(dict.lookupOrDefault<scalar>("rhoRef", 1)),
+    nSamples_(dict.lookupOrDefault<label>("N", 65536)),
+    fLower_(dict.lookupOrDefault<scalar>("fl", 25)),
+    fUpper_(dict.lookupOrDefault<scalar>("fu", 10000)),
+    startTime_(dict.lookupOrDefault<scalar>("startTime", 0)),
+    windowModelPtr_(windowModel::New(dict, nSamples_)),
+    graphFormat_(dict.lookupOrDefault<word>("graphFormat", "raw"))
 {
     // Check number of samples  - must be a power of 2 for our FFT
     bool powerOf2 = ((nSamples_ != 0) && !(nSamples_ & (nSamples_ - 1)));
diff --git a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H
index f7e739b7d0..28643cbf56 100644
--- a/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H
+++ b/src/randomProcesses/noise/noiseModels/noiseModel/noiseModel.H
@@ -30,7 +30,7 @@ Description
     Data is read from a dictionary, e.g.
 
     \verbatim
-    pRef            0;
+    rhoRef          0;
     N               4096;
     fl              25;
     fu              25;
@@ -40,10 +40,12 @@ Description
     where
     \table
         Property    | Description                   | Required  | Default value
-        pRef        | Reference pressure            | no        | 0
+        rhoRef      | Reference density             | no        | 1
         N           | Number of samples in sampling window | no | 65536 (2^16)
         fl          | Lower frequency bounds        | no        | 25
         fu          | Upper frequency bounds        | no        | 10000
+        startTime   | Start time                    | no        | 0
+        graphFormat | Graph format                  | no        | raw
     \endtable
 
 Note
@@ -94,8 +96,8 @@ protected:
         //- Copy of dictionary used for construction
         const dictionary dict_;
 
-        //- Reference pressure
-        scalar pRef_;
+        //- Reference density (to convert from kinematic to static pressure)
+        scalar rhoRef_;
 
         //- Number of samples in sampling window, default = 2^16
         label nSamples_;
@@ -111,7 +113,9 @@ protected:
 
         //- Window model
         autoPtr<windowModel> windowModelPtr_;
-        
+
+        //- Graph format
+        word graphFormat_;
 
 
     // Protected Member Functions
diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C
index 2014d721e3..fc03d03003 100644
--- a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C
+++ b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.C
@@ -43,7 +43,7 @@ addToRunTimeSelectionTable(noiseModel, pointNoise, dictionary);
 
 void pointNoise::filterTimeData
 (
-    const CSV<scalar>& pData,
+    const Function1Types::CSV<scalar>& pData,
     scalarField& t,
     scalarField& p
 )
@@ -80,12 +80,13 @@ void pointNoise::calculate()
 
     Info<< "Reading data file" << endl;
 
-    CSV<scalar> pData("pressure", dict_, "Data");
-    fileName baseFileName(pData.fName().lessExt());
+    Function1Types::CSV<scalar> pData("pressure", dict_, "Data");
 
     // Time and pressure history data
     scalarField t, p;
     filterTimeData(pData, t, p);
+    p *= rhoRef_;
+
     Info<< "    read " << t.size() << " values" << nl << endl;
 
     Info<< "Creating noise FFT" << endl;
@@ -93,38 +94,94 @@ void pointNoise::calculate()
 
     // Determine the windowing
     windowModelPtr_->validate(t.size());
+    const windowModel& win = windowModelPtr_();
+    const scalar deltaf = 1.0/(deltaT*win.nSamples());
+    fileName outDir(fileName("postProcessing")/"noise"/typeName);
 
     // Create the fft
     noiseFFT nfft(deltaT, p);
 
-    nfft -= pRef_;
 
-    graph Pf(nfft.RMSmeanPf(windowModelPtr_()));
-    Info<< "    Creating graph for " << Pf.title() << endl;
-    Pf.write(baseFileName + graph::wordify(Pf.title()), graphFormat_);
+    // Narrow band data
+    // ----------------
 
-    graph Lf(nfft.Lf(Pf));
-    Info<< "    Creating graph for " << Lf.title() << endl;
-    Lf.write(baseFileName + graph::wordify(Lf.title()), graphFormat_);
+    // RMS pressure [Pa]
+    graph Prmsf(nfft.RMSmeanPf(win));
+    Info<< "    Creating graph for " << Prmsf.title() << endl;
+    Prmsf.write(outDir, graph::wordify(Prmsf.title()), graphFormat_);
 
-    graph PSDf(nfft.PSDf(windowModelPtr_()));
+    // PSD [Pa^2/Hz]
+    graph PSDf(nfft.PSDf(win));
     Info<< "    Creating graph for " << PSDf.title() << endl;
-    PSDf.write(baseFileName + graph::wordify(PSDf.title()), graphFormat_);
+    PSDf.write(outDir, graph::wordify(PSDf.title()), graphFormat_);
 
-    graph PSD(nfft.PSD(PSDf));
-    Info<< "    Creating graph for " << PSD.title() << endl;
-    PSD.write(baseFileName + graph::wordify(PSD.title()), graphFormat_);
+    // PSD [dB/Hz]
+    graph PSDg
+    (
+        "PSD_dB_Hz(f)",
+        "f [Hz]",
+        "PSD(f) [dB_Hz]",
+        Prmsf.x(),
+        noiseFFT::PSD(PSDf.y())
+    );
+    Info<< "    Creating graph for " << PSDg.title() << endl;
+    PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_);
+
+    // SPL [dB]
+    graph SPLg
+    (
+        "SPL_dB(f)",
+        "f [Hz]",
+        "SPL(f) [dB]",
+        Prmsf.x(),
+        noiseFFT::SPL(PSDf.y()*deltaf)
+    );
+    Info<< "    Creating graph for " << SPLg.title() << endl;
+    SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
 
     labelList octave13BandIDs;
-    noiseFFT::octaveFrequenciesIDs(Pf.x(), fLower_, fUpper_, 3, octave13BandIDs);
+    scalarField octave13FreqCentre;
+    noiseFFT::octaveBandInfo
+    (
+        Prmsf.x(),
+        fLower_,
+        fUpper_,
+        3,
+        octave13BandIDs,
+        octave13FreqCentre
+    );
 
-    graph Ldelta(nfft.Ldelta(Lf, octave13BandIDs));
-    Info<< "    Creating graph for " << Ldelta.title() << endl;
-    Ldelta.write(baseFileName + graph::wordify(Ldelta.title()), graphFormat_);
 
-    graph Pdelta(nfft.Pdelta(Pf, octave13BandIDs));
-    Info<< "    Creating graph for " << Pdelta.title() << endl;
-    Pdelta.write(baseFileName + graph::wordify(Pdelta.title()), graphFormat_);
+    // 1/3 octave data
+    // ---------------
+
+    // PSD [Pa^2/Hz]
+    graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false));
+
+    // Integrated PSD = P(rms)^2 [Pa^2]
+    graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true));
+
+    graph PSD13g
+    (
+        "PSD13_dB_Hz(fm)",
+        "fm [Hz]",
+        "PSD(fm) [dB_Hz]",
+        octave13FreqCentre,
+        noiseFFT::PSD(PSD13f.y())
+    );
+    Info<< "    Creating graph for " << PSD13g.title() << endl;
+    PSD13g.write(outDir, graph::wordify(PSD13g.title()), graphFormat_);
+
+    graph SPL13g
+    (
+        "SPL13_dB(fm)",
+        "fm [Hz]",
+        "SPL(fm) [dB]",
+        octave13FreqCentre,
+        noiseFFT::SPL(Prms13f2.y())
+    );
+    Info<< "    Creating graph for " << SPL13g.title() << endl;
+    SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_);
 }
 
 
@@ -132,8 +189,7 @@ void pointNoise::calculate()
 
 pointNoise::pointNoise(const dictionary& dict)
 :
-    noiseModel(dict),
-    graphFormat_(dict.lookupOrDefault<word>("graphFormat", "raw"))
+    noiseModel(dict)
 {}
 
 
diff --git a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H
index 020e462238..30bc669ae9 100644
--- a/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H
+++ b/src/randomProcesses/noise/noiseModels/pointNoise/pointNoise.H
@@ -37,10 +37,10 @@ Description
     // Must be a power of 2, default = 2^16 (=65536)
     N               4096;
 
-    // Lower frequency bounds 
+    // Lower frequency bounds
     fl              25;
 
-    // Upper frequency bounds 
+    // Upper frequency bounds
     fu              25;
 
     // Start time
@@ -100,17 +100,11 @@ class pointNoise
 
 protected:
 
-    // Protected data
-
-        //- Graph format
-        word graphFormat_;
-
-
     // Protected Member Functions
 
         void filterTimeData
         (
-            const CSV<scalar>& pData,
+            const Function1Types::CSV<scalar>& pData,
             scalarField& t,
             scalarField& p
         );
diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
index 1dcc319766..632fda9d72 100644
--- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
+++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C
@@ -46,6 +46,8 @@ addToRunTimeSelectionTable(noiseModel, surfaceNoise, dictionary);
 
 void surfaceNoise::initialise(const dictionary& dict)
 {
+    label nAvailableTimes = 0;
+
     // All reading performed on the master processor only
     if (Pstream::master())
     {
@@ -81,10 +83,22 @@ void surfaceNoise::initialise(const dictionary& dict)
         startTimeIndex_ = findStartTimeIndex(allTimes, startTime_);
 
         // Determine the windowing
-        label nAvailableTimes = allTimes.size() - startTimeIndex_;
-        label nRequiredTimes = windowModelPtr_->validate(nAvailableTimes);
+        nAvailableTimes = allTimes.size() - startTimeIndex_;
+    }
 
+    Pstream::scatter(pIndex_);
+    Pstream::scatter(startTimeIndex_);
+    Pstream::scatter(nAvailableTimes);
+
+
+    // Note: all processors should call the windowing validate function
+    label nRequiredTimes = windowModelPtr_->validate(nAvailableTimes);
+
+    if (Pstream::master())
+    {
         // Restrict times
+        const instantList allTimes = readerPtr_->times();
+
         times_.setSize(nRequiredTimes);
         forAll(times_, timeI)
         {
@@ -92,13 +106,12 @@ void surfaceNoise::initialise(const dictionary& dict)
         }
         deltaT_ = checkUniformTimeStep(times_);
 
+        // Read the surface geometry
         const meshedSurface& surf = readerPtr_->geometry();
         nFace_ = surf.size();
     }
 
-    Pstream::scatter(pIndex_);
     Pstream::scatter(times_);
-    Pstream::scatter(startTimeIndex_);
     Pstream::scatter(deltaT_);
     Pstream::scatter(nFace_);
 }
@@ -154,6 +167,9 @@ void surfaceNoise::readSurfaceData
                 // Read pressure at all faces for time timeI
                 scalarField p(readerPtr_->field(timeI, pIndex_, scalar(0)));
 
+                // Apply conversions
+                p *= rhoRef_;
+
                 // Send subset of faces to each processor
                 for (label procI = 0; procI < Pstream::nProcs(); procI++)
                 {
@@ -201,7 +217,7 @@ void surfaceNoise::readSurfaceData
 
             forAll(p, faceI)
             {
-                pData[faceI][i] = p[faceI];
+                pData[faceI][i] = p[faceI]*rhoRef_;
             }
         }
     }
@@ -214,7 +230,7 @@ void surfaceNoise::readSurfaceData
 }
 
 
-void surfaceNoise::writeSurfaceData
+Foam::scalar surfaceNoise::writeSurfaceData
 (
     const word& fName,
     const word& groupName,
@@ -245,6 +261,7 @@ void surfaceNoise::writeSurfaceData
 
         pBufs.finishedSends();
 
+        scalar areaAverage = 0;
         if (Pstream::master())
         {
             const meshedSurface& surf = readerPtr_->geometry();
@@ -278,7 +295,20 @@ void surfaceNoise::writeSurfaceData
                 allData,
                 false
             );
+
+            // TODO: Move faceAreas to demand-driven function in MeshedSurface
+            scalarField faceAreas(surf.faces().size());
+            forAll(faceAreas, i)
+            {
+                faceAreas[i] = surf.faces()[i].mag(surf.points());
+            }
+
+//            areaAverage = sum(allData*faceAreas)/sum(faceAreas);
+            areaAverage = sum(allData)/allData.size();
         }
+        Pstream::scatter(areaAverage);
+
+        return areaAverage;
     }
     else
     {
@@ -294,6 +324,91 @@ void surfaceNoise::writeSurfaceData
             data,
             false
         );
+
+        // TODO: Move faceAreas to demand-driven function in MeshedSurface
+        scalarField faceAreas(surf.faces().size());
+        forAll(faceAreas, i)
+        {
+            faceAreas[i] = surf.faces()[i].mag(surf.points());
+        }
+
+//        return sum(data*faceAreas)/sum(faceAreas);
+        return sum(data)/data.size();
+    }
+}
+
+
+Foam::scalar surfaceNoise::surfaceAverage
+(
+    const scalarField& data,
+    const labelList& procFaceOffset
+) const
+{
+    if (Pstream::parRun())
+    {
+        // Collect the surface data so that we can output the surfaces
+
+        PstreamBuffers pBufs(Pstream::nonBlocking);
+
+        if (!Pstream::master())
+        {
+            UOPstream toProc(0, pBufs);
+            toProc << data;
+        }
+
+        pBufs.finishedSends();
+
+        scalar areaAverage = 0;
+        if (Pstream::master())
+        {
+            const meshedSurface& surf = readerPtr_->geometry();
+
+            scalarField allData(surf.size());
+
+            forAll(data, faceI)
+            {
+                // Master procFaceOffset is zero...
+                allData[faceI] = data[faceI];
+            }
+
+            for (label procI = 1; procI < Pstream::nProcs(); procI++)
+            {
+                UIPstream fromProc(procI, pBufs);
+                scalarList dataSlice(fromProc);
+                forAll(dataSlice, i)
+                {
+                    label faceI = procFaceOffset[procI] + i;
+                    allData[faceI] = dataSlice[i];
+                }
+            }
+
+            // TODO: Move faceAreas to demand-driven function in MeshedSurface
+            scalarField faceAreas(surf.faces().size());
+            forAll(faceAreas, i)
+            {
+                faceAreas[i] = surf.faces()[i].mag(surf.points());
+            }
+
+//            areaAverage = sum(allData*faceAreas)/sum(faceAreas);
+            areaAverage = sum(allData)/allData.size();
+        }
+        Pstream::scatter(areaAverage);
+
+        return areaAverage;
+    }
+    else
+    {
+        const meshedSurface& surf = readerPtr_->geometry();
+
+        // TODO: Move faceAreas to demand-driven function in MeshedSurface
+        scalarField faceAreas(surf.faces().size());
+        forAll(faceAreas, i)
+        {
+            faceAreas[i] = surf.faces()[i].mag(surf.points());
+        }
+
+//        return sum(data*faceAreas)/sum(faceAreas);
+        return sum(data)/data.size();
     }
 }
 
@@ -358,123 +473,226 @@ void surfaceNoise::calculate()
     const label nLocalFace = pData.size();
     const scalarField freq1(noiseFFT::frequencies(nSamples_, deltaT_));
     const label nFFT = freq1.size()/fftWriteInterval_;
-    List<scalarField> surfPf(nFFT);
-    List<scalarField> surfLf(nFFT);
-    List<scalarField> surfPSD(nFFT);
-    forAll(surfPf, freqI)
+    List<scalarField> surfPrmsf(nFFT);
+    List<scalarField> surfPSDf(nFFT);
+    forAll(surfPrmsf, freqI)
     {
-        surfPf[freqI].setSize(nLocalFace);
-        surfLf[freqI].setSize(nLocalFace);
-        surfPSD[freqI].setSize(nLocalFace);
+        surfPrmsf[freqI].setSize(nLocalFace);
+        surfPSDf[freqI].setSize(nLocalFace);
     }
 
     // Storage for 1/3 octave data
     labelList octave13BandIDs;
-    noiseFFT::octaveFrequenciesIDs(freq1, fLower_, fUpper_, 3, octave13BandIDs);
+    scalarField octave13FreqCentre;
+    noiseFFT::octaveBandInfo
+    (
+        freq1,
+        fLower_,
+        fUpper_,
+        3,
+        octave13BandIDs,
+        octave13FreqCentre
+    );
 
-    List<scalarField> surfPdelta(octave13BandIDs.size() - 1);
-    List<scalarField> surfLdelta(octave13BandIDs.size() - 1);
-    forAll(surfPdelta, freqI)
+    List<scalarField> surfPSD13f(octave13BandIDs.size() - 1);
+    List<scalarField> surfPrms13f2(octave13BandIDs.size() - 1);
+    forAll(surfPSD13f, freqI)
     {
-        surfPdelta[freqI].setSize(nLocalFace);
-        surfLdelta[freqI].setSize(nLocalFace);
+        surfPSD13f[freqI].setSize(nLocalFace);
+        surfPrms13f2[freqI].setSize(nLocalFace);
     }
 
+    const windowModel& win = windowModelPtr_();
+
     forAll(pData, faceI)
     {
         const scalarField& p = pData[faceI];
 
         noiseFFT nfft(deltaT_, p);
-
-        nfft -= pRef_;
-        graph Pf(nfft.RMSmeanPf(windowModelPtr_()));
-        graph Lf(nfft.Lf(Pf));
-        graph PSDf(nfft.PSDf(windowModelPtr_()));
-        graph PSD(nfft.PSD(PSDf));
+        graph Prmsf(nfft.RMSmeanPf(win));
+        graph PSDf(nfft.PSDf(win));
 
         // Store the frequency results in slot for face of surface
-        forAll(surfPf, i)
+        forAll(surfPrmsf, i)
         {
             label freqI = (i + 1)*fftWriteInterval_ - 1;
-            surfPf[i][faceI] = Pf.y()[freqI];
-            surfLf[i][faceI] = Lf.y()[freqI];
-            surfPSD[i][faceI] = PSD.y()[freqI];
+            surfPrmsf[i][faceI] = Prmsf.y()[freqI];
+            surfPSDf[i][faceI] = PSDf.y()[freqI];
         }
 
-        graph Pdelta(nfft.Pdelta(Pf, octave13BandIDs));
-        graph Ldelta(nfft.Ldelta(Lf, octave13BandIDs));
+        // PSD [Pa^2/Hz]
+        graph PSD13f(nfft.octaves(PSDf, octave13BandIDs, false));
+
+        // Integrated PSD = P(rms)^2 [Pa^2]
+        graph Prms13f2(nfft.octaves(PSDf, octave13BandIDs, true));
 
         // Store the 1/3 octave results in slot for face of surface
-        forAll(surfPdelta, freqI)
+        forAll(surfPSD13f, freqI)
         {
-            surfPdelta[freqI][faceI] = Pdelta.y()[freqI];
-            surfLdelta[freqI][faceI] = Ldelta.y()[freqI];
+            surfPSD13f[freqI][faceI] = PSD13f.y()[freqI];
+            surfPrms13f2[freqI][faceI] = Prms13f2.y()[freqI];
         }
 
         // Free the storage for p
 //        p.clear();
     }
 
-    Info<< "Writing fft surface data" << endl;
+    // Output directory for graphs
+    fileName outDir(fileName("postProcessing")/"noise"/typeName);
 
-    forAll(surfPf, i)
+    const scalar deltaf = 1.0/(deltaT_*win.nSamples());
+    Info<< "Writing fft surface data" << endl;
     {
-        label freqI = i*fftWriteInterval_;
-        const word& fName = inputFileName_.name(true);
-        const word gName = "fft";
-        writeSurfaceData
+        scalarField PrmsfAve(surfPrmsf.size(), 0);
+        scalarField PSDfAve(surfPrmsf.size(), 0);
+        scalarField fOut(surfPrmsf.size(), 0);
+
+        forAll(surfPrmsf, i)
+        {
+            label freqI = i*fftWriteInterval_;
+            fOut[i] = freq1[freqI];
+            const word& fName = inputFileName_.name(true);
+            const word gName = "fft";
+            PrmsfAve[i] = writeSurfaceData
+            (
+                fName,
+                gName,
+                "Prmsf",
+                freq1[freqI],
+                surfPrmsf[i],
+                procFaceOffset
+            );
+
+            PSDfAve[i] = writeSurfaceData
+            (
+                fName,
+                gName,
+                "PSDf",
+                freq1[freqI],
+                surfPSDf[i],
+                procFaceOffset
+            );
+            writeSurfaceData
+            (
+                fName,
+                gName,
+                "PSD",
+                freq1[freqI],
+                noiseFFT::PSD(surfPSDf[i]),
+                procFaceOffset
+            );
+            writeSurfaceData
+            (
+                fName,
+                gName,
+                "SPL",
+                freq1[freqI],
+                noiseFFT::SPL(surfPSDf[i]*deltaf),
+                procFaceOffset
+            );
+        }
+
+        graph Prmsfg
         (
-            fName,
-            gName,
-            "Pf",
-            freq1[freqI],
-            surfPf[i],
-            procFaceOffset
+            "Average Prms(f)",
+            "f [Hz]",
+            "P(f) [Pa]",
+            fOut,
+            PrmsfAve
         );
-        writeSurfaceData
+        Prmsfg.write(outDir, graph::wordify(Prmsfg.title()), graphFormat_);
+
+        graph PSDfg
         (
-            fName,
-            gName,
-            "Lf",
-            freq1[freqI],
-            surfLf[i],
-            procFaceOffset
+            "Average PSD_f(f)",
+            "f [Hz]",
+            "PSD(f) [PaPa_Hz]",
+            fOut,
+            PSDfAve
         );
-        writeSurfaceData
+        PSDfg.write(outDir, graph::wordify(PSDfg.title()), graphFormat_);
+
+        graph PSDg
         (
-            fName,
-            gName,
-            "PSD",
-            freq1[freqI],
-            surfPSD[i],
-            procFaceOffset
+            "Average PSD_dB_Hz(f)",
+            "f [Hz]",
+            "PSD(f) [dB_Hz]",
+            fOut,
+            noiseFFT::PSD(PSDfAve)
         );
+        PSDg.write(outDir, graph::wordify(PSDg.title()), graphFormat_);
+
+        graph SPLg
+        (
+            "Average SPL_dB(f)",
+            "f [Hz]",
+            "SPL(f) [dB]",
+            fOut,
+            noiseFFT::SPL(PSDfAve*deltaf)
+        );
+        SPLg.write(outDir, graph::wordify(SPLg.title()), graphFormat_);
     }
 
-    Info<< "Writing one-third octave surface data" << endl;
 
-    forAll(surfPdelta, i)
+    Info<< "Writing one-third octave surface data" << endl;
     {
-        const word& fName = inputFileName_.name(true);
-        const word gName = "oneThirdOctave";
-        writeSurfaceData
+        scalarField PSDfAve(surfPSD13f.size(), 0);
+        scalarField Prms13f2Ave(surfPSD13f.size(), 0);
+
+        forAll(surfPSD13f, i)
+        {
+            const word& fName = inputFileName_.name(true);
+            const word gName = "oneThirdOctave";
+            PSDfAve[i] = writeSurfaceData
+            (
+                fName,
+                gName,
+                "PSD13f",
+                octave13FreqCentre[i],
+                surfPSD13f[i],
+                procFaceOffset
+            );
+            writeSurfaceData
+            (
+                fName,
+                gName,
+                "PSD13",
+                octave13FreqCentre[i],
+                noiseFFT::PSD(surfPSD13f[i]),
+                procFaceOffset
+            );
+            writeSurfaceData
+            (
+                fName,
+                gName,
+                "SPL13",
+                octave13FreqCentre[i],
+                noiseFFT::SPL(surfPrms13f2[i]),
+                procFaceOffset
+            );
+
+            Prms13f2Ave[i] = surfaceAverage(surfPrms13f2[i], procFaceOffset);
+        }
+
+        graph PSD13g
         (
-            fName,
-            gName,
-            "Pdelta",
-            octave13BandIDs[i],
-            surfPdelta[i],
-            procFaceOffset
+            "Average PSD13_dB_Hz(fm)",
+            "fm [Hz]",
+            "PSD(fm) [dB_Hz]",
+            octave13FreqCentre,
+            noiseFFT::PSD(PSDfAve)
         );
-        writeSurfaceData
+        PSD13g.write(outDir, graph::wordify(PSD13g.title()), graphFormat_);
+
+        graph SPL13g
         (
-            fName,
-            gName,
-            "Ldelta",
-            octave13BandIDs[i],
-            surfLdelta[i],
-            procFaceOffset
+            "Average SPL13_dB(fm)",
+            "fm [Hz]",
+            "SPL(fm) [dB]",
+            octave13FreqCentre,
+            noiseFFT::SPL(Prms13f2Ave)
         );
+        SPL13g.write(outDir, graph::wordify(SPL13g.title()), graphFormat_);
     }
 }
 
@@ -485,6 +703,3 @@ void surfaceNoise::calculate()
 } // End namespace Foam
 
 // ************************************************************************* //
-
-
-
diff --git a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H
index e2416e838a..3f80ee8b4a 100644
--- a/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H
+++ b/src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.H
@@ -158,7 +158,8 @@ protected:
         );
 
         //- Write surface data to file
-        void writeSurfaceData
+        //  Returns the area average value
+        scalar writeSurfaceData
         (
             const word& fName,
             const word& groupName,
@@ -168,6 +169,12 @@ protected:
             const labelList& procFaceOffset
         ) const;
 
+        //- Calculate the area average value
+        scalar surfaceAverage
+        (
+            const scalarField& data,
+            const labelList& procFaceOffset
+        ) const;
 
 
 public: