Having issue with units for wav file spectral density analysis (db artificially low)

Question

I have the below code which is adapted from other code I found here. It runs fine and the plot looks as expected except the dB levels are about 80dB too low. Note, the intended unit is dB SPL. I would appreciate any input on fixing this, I am thinking it has something to do with needing to include the 1uPa reference but I can't track back the issue and I'm just a poor biologist trying to learn to code!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 96000;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);
% spectrogram dB scale
spectrogram_dB_scale = 100;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
gain = 20
sensitivity = - 165
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[data,Fs]=audioread('sample.wav'); %(newer matlab)
samples = length(data);
dt = 1/Fs;
t = (0:dt:(samples-1)*dt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(data,w,NOVERLAP,NFFT,Fs);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum) - sensitivity - gain ;


figure(1),semilogx(freq,sensor_spectrum_dB);grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB (L))');

Answer 1

1. Regardless, you want to plot the output of pwelch using 10*log10 since you're dealing with a power spectrum.

2. Without knowing a little more about your data it's going to be hard to help... is it sound data? physiological? what is the sensitivity unit? What is gain ??

3. Scaling the PSD