Average of multiple Thompson Multitaper Power Spectral Density Estimates

I am doing signal processing of fricative sounds. I have managed to plot the tokens I have on the same plot, but I also want to plot the average spectra, based on the other PSDs . Of course, I would like it if I could also plot the average spectra in a different color! Any help would be appreciated.

The current spectra I have were created using the functions:

Hs = spectrum.mtm;
fs = 44100 %this is the recording quality of the signal;
x = wavread('filename');
y = wavread('otherfilename');
z = wavread('lastfilename');
x1 = psd(Hs, x, 'fs', fs);
y2 = psd(Hs, y, 'fs', fs);
z1 = psd(Hs, z, 'fs', fs);

hold on;
plot(x1);
plot(y1);
plot(z1);

I want to plot the average spectrum on top of those spectrums (in black, if possible)! Again, thank you for any help!

I have attached the .wav files in a .zip folder.

Phil

[Image Analyst]

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

Answer 1

[Image Analyst]

So can't you just take the mean and plot it? I mean

meanSpectrum = (x1+y1+z1)/3;
plot(meanSpectrum, 'k-', 'LineWidth', 3);
legend('x1', 'y1', 'z1', 'meanSpectrum');

It seems really trivial, or am I missing something???

[Phil]

I'm a total newbie to MATLAB, so simple things often take a while of searching through help forums.

Anyway, I got an error from the first line of code, it said:

??? Undefined function or method 'plus' for input arguments of type 'dspdata.psd'.

What I just tried instead was the following:

    meanSpectrum = (x, y, z)/3;
    mean = (Hs, meanSpectrum, 'fs', fs);
    plot(mean)

would this give me an average as well?

I do have one issue though, when I try to add the specifications for line width and color, I get the following error:

>> plot(mean, 'k-', 'LineWidth', 3) ??? Error using ==> dspdata.abstractfreqresp.plot Too many input arguments.

Thanks again for the initial help and any further help!

[Image Analyst]

Well almost, except that you called your variables x1, y1, and z1 instead of x, y, and z. Plus you can't use (or should NOT use) mean as the name of one of your variables because it's the name of a built in function. I don't know the type of x - I've not used that function. What does this show

whos x
x

[Phil]

x y z are the names of the raw wave files I loaded into MATLAB before the I did the multitaper spectral density estimate. Basically I averaged all the waves and then ran a multitaper spectral density estimate on the average wave. I'm not sure this works well, because I think the average then becomes the lowest spectrum , when it should be somewhere in the center of the other spectrums. I think what I want to do is average the multitaper spectrums and then plot that mean spectrum.

>> whos x

Name         Size            Bytes  Class     Attributes

    x         1323x1             10584  double

if I punch in x it gives me the vectors of the wave I loaded in.

[Image Analyst]

But x, y, and z are not the spectrum. They're in the time domain. Your spectra are the badly-named x1, y1, and z1. Are you sure the length of x, y, and z are all the same and the lengths of x1, y1, and z1 are all the same? Probably because you say you were able to average the wave signals. But since the Fourier transform is a linear transform the transform of the sum is the sum of the transforms so you should get something in the middle, not just exactly the lowest spectrum. Why don't you upload your wave files?

[Phil]

I attached a .zip file of the .wav files I am using.

The issue is that once I do a Thompson multitaper power spectral density estimate, it won't let me add the .psd in order to average. Would an ensemble average work? How would I implement it if it would?

Thanks for the help so far.

[Image Analyst]

Try this:

Hs = spectrum.mtm;
fs = 44100 %this is the recording quality of the signal;
x = wavread('C:\Users\Phil\Documents\sheppard_1.wav');
y = wavread('C:\Users\Phil\Documents\sheppard_2.wav');
z = wavread('C:\Users\Phil\Documents\sheppard_3.wav');
xPsd = psd(Hs, x, 'fs', fs);
yPsd = psd(Hs, y, 'fs', fs);
zPsd = psd(Hs, z, 'fs', fs);

hold on;
plot(xPsd.Data, 'r-', 'LineWidth', 2);
plot(yPsd.Data, 'g-', 'LineWidth', 2);
plot(zPsd.Data, 'b-', 'LineWidth', 2);
meanPsd = (xPsd.Data + yPsd.Data + zPsd.Data)/3;
plot(meanPsd, 'k-', 'LineWidth', 2);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
grid on;
legend('x', 'y', 'z', 'Mean');

[Phil]

This does indeed get an average, but the plot is much different than what I had in mind. I am trying to get the spectra on a plot where the x-axis is frequency and the y-axis is power/frequency (dB/Hz)

I attached an example plot in .bmp and .fig format. The only thing I am missing is the mean. I'm content with everything being the same color, except I want to plot the mean PSD in bold black.

[Image Analyst]

xPSD has another structure member that is the frequencies. So just plot it:

plot(xPsd.Frequencies, meanPsd, 'k-', 'LineWidth', 2);

If this answered your question, please go ahead and officially mark the answer as "Accepted" - Thanks.

[Phil]

Unfortunately, this doesn't appear to work. What I get if I enter the data as:

meanPsd = (xPsd.Data + yPsd.Data + zPsd.Data)/3;

is a flat line across the y-axis at 0. If I enter the meanPsd as:

meanPsd = (xPsd.Frequencies + yPsd.Frequencies + zPsd.Frequencies)/3;

I get a diagonal line from 0 to just over 20. What I am trying to get is a mean that resembles the spectrums, I attached an image from another source. It has some extra stuff I can worry about later, but the mean for that particular set of data is a black line that resembles the other spractra.

[Image Analyst]

Who told you to average the frequencies to get the PSD? Not me. I wouldn't do that since that's nonsense. Look at what I told you. Alright, see this code:

Hs = spectrum.mtm;
fs = 44100 %this is the recording quality of the signal;
x = wavread('C:\Users\Mark\Documents\Temporary\sheppard_1.wav');
y = wavread('C:\Users\Mark\Documents\Temporary\sheppard_2.wav');
z = wavread('C:\Users\Mark\Documents\Temporary\sheppard_3.wav');
xPsd = psd(Hs, x, 'fs', fs);
yPsd = psd(Hs, y, 'fs', fs);
zPsd = psd(Hs, z, 'fs', fs);

hold on;
plot(xPsd.Frequencies, xPsd.Data, 'r-', 'LineWidth', 2);
plot(yPsd.Frequencies, yPsd.Data, 'g-', 'LineWidth', 2);
plot(zPsd.Frequencies, zPsd.Data, 'b-', 'LineWidth', 2);
meanPsd = (xPsd.Data + yPsd.Data + zPsd.Data)/3;
plot(xPsd.Frequencies, meanPsd, 'k-', 'LineWidth', 5);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
grid on;
legend('x', 'y', 'z', 'Mean');

[Phil]

I'm clearly not getting something you told me to do and I tried that frequencies in order to solve the issue, but it didn't help. Here is what I have, because when I do what you suggest, I get a flat line at 0 on the y-axis:

a = wavread('sheppard_1');
b = wavread('sheppard_2');
c = wavread('sheppard_3');
d = wavread('sheppard_4');
e = wavread('sheppard_5');
f = wavread('sheppard_6');
g = wavread('sheppard_7');
h = wavread('sheppard_8');
i = wavread('sheppard_9');

Hs = spectrum.mtm;
fs = 44100;

apsd = psd(Hs, a, 'fs', fs');
bpsd = psd(Hs, b, 'fs', fs');
cpsd = psd(Hs, c, 'fs', fs');
dpsd = psd(Hs, d, 'fs', fs');
epsd = psd(Hs, e, 'fs', fs');
fpsd = psd(Hs, f, 'fs', fs');
gpsd = psd(Hs, g, 'fs', fs');
hpsd = psd(Hs, h, 'fs', fs');
ipsd = psd(Hs, i, 'fs', fs');

meanpsd = (apsd.data + bpsd.data + cpsd.data + dpsd.data + epsd.data + fpsd.data + gpsd.data + hpsd.data + ipsd.data)/9;

hold on

plot(apsd);
plot(bpsd);
plot(cpsd);
plot(dpsd);
plot(epsd);
plot(fpsd);
plot(gpsd);
plot(hpsd);
plot(ipsd);
plot(apsd.Frequencies, meanpsd, 'k-', 'LineWidth', 2);

hold off

I attached an image of the results.

[Phil]

Right, I see what you've done, but that's not what I'm trying to do. I'm trying to create a plot like the one I just attached, but with a mean as well.

[Image Analyst]

Mine does that - the main difference is that I used different colors so that you can see the difference between the different curves, and I used black for the mean (like you wanted but didn't do yourself) and I plotted the mean in a thicker line so it is more noticeable.

Plus you clearly have code you haven't shared. For example there is code where you set the x and y labels. I can see it in your results plot but not in your code so some code is missing. Plus the y scale of your plots are different than mine.

[Phil]

I shared the exact code that I used to create the plot that I attached. I see your plot, but it isn't shaped like a "spectral slice" like mine is, which is what I'm trying to do!

Ok, I tried the following and it is giving me what I want, but the x-axis is no longer Frequency(kHz):

a = wavread('sheppard_1');
b = wavread('sheppard_2');
c = wavread('sheppard_3');
d = wavread('sheppard_4');
e = wavread('sheppard_5');
f = wavread('sheppard_6');
g = wavread('sheppard_7');
h = wavread('sheppard_8');
i = wavread('sheppard_9');

fs = 44100;

[pxxa,fa] = pmtm(a, 3, length(a), fs);
[pxxb,fb] = pmtm(b, 3, length(b), fs);
[pxxc,fc] = pmtm(c, 3, length(c), fs);
[pxxd,fd] = pmtm(d, 3, length(d), fs);
[pxxe,fe] = pmtm(e, 3, length(e), fs);
[pxxf,ff] = pmtm(f, 3, length(f), fs);
[pxxg,fg] = pmtm(g, 3, length(g), fs);
[pxxh,fh] = pmtm(h, 3, length(h), fs);
[pxxi,fi] = pmtm(i, 3, length(i), fs);

meanpxx = (pxxa + pxxb + pxxc + pxxd + pxxe + pxxf + pxxg + pxxh + pxxi)/9;

hold on

plot(fa, 10*log10(pxxa));
plot(fb, 10*log10(pxxb));
plot(fc, 10*log10(pxxc));
plot(fd, 10*log10(pxxd));
plot(fe, 10*log10(pxxe));
plot(ff, 10*log10(pxxf));
plot(fg, 10*log10(pxxg));
plot(fh, 10*log10(pxxh));
plot(fi, 10*log10(pxxi));
plot(fa, 10*log10(meanpxx), 'k-', 'LineWidth', 2);

hold off

I attached the result, the only issue is the scale is a bit weird, it reads as x-axis value X 10^4, which is not what I want...

[Image Analyst]

This is my code. Originally you were plotting the first 9 as a function of element number, not frequency, and only the mean PSD was vs. frequency, so it was all messed up.

Hs = spectrum.mtm;
fs = 44100 %this is the recording quality of the signal;
w1 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_1.wav');
w2 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_2.wav');
w3 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_3.wav');
w4 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_4.wav');
w5 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_5.wav');
w6 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_6.wav');
w7 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_7.wav');
w8 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_8.wav');
w9 = wavread('C:\Users\Mark\Documents\Temporary\sheppard_9.wav');
psd1 = psd(Hs, w1, 'fs', fs);
psd2 = psd(Hs, w2, 'fs', fs);
psd3 = psd(Hs, w3, 'fs', fs);
psd4 = psd(Hs, w4, 'fs', fs);
psd5 = psd(Hs, w5, 'fs', fs);
psd6 = psd(Hs, w6, 'fs', fs);
psd7 = psd(Hs, w7, 'fs', fs);
psd8 = psd(Hs, w8, 'fs', fs);
psd9 = psd(Hs, w9, 'fs', fs);

hold on;
plot(psd1.Frequencies, psd1.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd2.Frequencies, psd2.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd3.Frequencies, psd3.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd4.Frequencies, psd4.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd5.Frequencies, psd5.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd6.Frequencies, psd6.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd7.Frequencies, psd7.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd8.Frequencies, psd8.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
plot(psd9.Frequencies, psd9.Data, '-', 'Color', rand(1,3), 'LineWidth', 2);
meanPsd = (psd1.Data + psd2.Data + psd3.Data + psd4.Data + psd5.Data + psd6.Data + psd7.Data + psd8.Data + psd9.Data)/9;
plot(psd1.Frequencies, meanPsd, 'k-', 'LineWidth', 5);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
grid on;
legend('x', 'y', 'z', 'Mean');

[Phil]

Based on what you were doing I changed the functions to:

[pxxa,fa] = pmtm(a, 3, length(a), fs);
[pxxb,fb] = pmtm(b, 3, length(b), fs);
[pxxc,fc] = pmtm(c, 3, length(c), fs);
[pxxd,fd] = pmtm(d, 3, length(d), fs);
[pxxe,fe] = pmtm(e, 3, length(e), fs);
[pxxf,ff] = pmtm(f, 3, length(f), fs);
[pxxg,fg] = pmtm(g, 3, length(g), fs);
[pxxh,fh] = pmtm(h, 3, length(h), fs);
[pxxi,fi] = pmtm(i, 3, length(i), fs);

Instead of the psd function! I'm not sure why it worked, but it did, and then I was able to follow your original suggestion of just adding them and dividing by the total, which gave me a plot like the one I attached, which is exactly what I'm trying to plot. To get the axis the way I wanted, I had to manually change it. I'm gonna try and figure out some code to automate that as well.

Thanks for taking the time to help me! I'm sure I couldn't have figured this out without your suggestions.

Thanks!