Smoothing Power Spectral Density

Matlab. Having influenced by signal smoothing algorithms like moving average algorithm, I experimentally developed a smoothing algorithm for total of power spectral density. My purpose was to smooth delusive sharp variations in the graph so that I can do some calculations better. It balances between sharp peaks and falls. My first aim is to see the attack points better by this algorithm.

It takes the total of power spectral density to keep all frequencies in the instant time. First step takes average value of the previous and next points. Then, it adds this value as current and next values. Second step does the same calculation, however it takes the average value of the previous and next two jump points of first smoothed graph. Both steps have the a length correction calculation. I used two piano segments for the demonstration: piano1 and piano2. Please see the graphs for two segments below.

function psd_smoothing
tic
[y, Fs] = audioread('piano.wav');         %input audio file (mono)

[~,~,T,P]=spectrogram(y,256,[],[],Fs);

totalpower = sum(P);

%PSD GRAPH SMOOTHING
 
%Applying First Step
smooth1 = totalpower(1);
for i = 2:2:length(T)-2
   new1 = (totalpower(i-1)+totalpower(i+1))/2;
   smooth1 = [smooth1, new1, new1]; 
end

if abs(length(T) - length(smooth1)) > 0
    for j = abs(length(T) - length(smooth1)):-1:1
        addition = length(T)+1-j;
        smooth1 = [smooth1, totalpower(addition)];
    end
end    

%Applying Second Step
smooth2 = smooth1(1);
for i = 2:2:length(T)-2
   new2 = (smooth1(i-1)+smooth1(i+2))/2;
   smooth2 = [smooth2, new2, new2]; 
end

if abs(length(T) - length(smooth2)) > 0
    for j = abs(length(T) - length(smooth2)):-1:1
        addition = length(T)+1-j;
        smooth2 = [smooth2, smooth1(addition)];
    end
end

%Plotting Stage
subplot(3,1,1);
plot(T,P);
subplot(3,1,2);
plot(T,totalpower);
subplot(3,1,3);
plot(T,smooth2);

toc
end   %end of main function

function psd_smoothing

tic

[y, Fs] = audioread('piano.wav'); %input audio file (mono)

[~,~,T,P]=spectrogram(y,256,[],[],Fs);

totalpower = sum(P);

%PSD GRAPH SMOOTHING

%Applying First Step

smooth1 = totalpower(1);

for i = 2:2:length(T)-2

new1 = (totalpower(i-1)+totalpower(i+1))/2;

smooth1 = [smooth1, new1, new1];

end

if abs(length(T) - length(smooth1)) > 0

for j = abs(length(T) - length(smooth1)):-1:1

addition = length(T)+1-j;

smooth1 = [smooth1, totalpower(addition)];

end

%Applying Second Step

smooth2 = smooth1(1);

for i = 2:2:length(T)-2

new2 = (smooth1(i-1)+smooth1(i+2))/2;

smooth2 = [smooth2, new2, new2];

end

if abs(length(T) - length(smooth2)) > 0

for j = abs(length(T) - length(smooth2)):-1:1

addition = length(T)+1-j;

smooth2 = [smooth2, smooth1(addition)];

end

%Plotting Stage

subplot(3,1,1);

plot(T,P);

subplot(3,1,2);

plot(T,totalpower);

subplot(3,1,3);

plot(T,smooth2);

toc

end %end of main function

References:

[1] Power Spectral Density, Matlab Documentation

[2] Signal Smoothing Algorithms