Average of multiple Thompson Multitaper Power Spectral Density Estimates

Question

Phil 2014-8-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/144481-average-of-multiple-thompson-multitaper-power-spectral-density-estimates

编辑： Image Analyst 2014-8-3

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Image Analyst 2014-8-2

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Image Analyst 2014-8-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/144481-average-of-multiple-thompson-multitaper-power-spectral-density-estimates#answer_147212

在 MATLAB Online 中打开

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???

16 个评论
显示 14更早的评论隐藏 14更早的评论

Phil 2014-8-3

在 MATLAB Online 中打开

example.bmp

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 2014-8-3

在 MATLAB Online 中打开

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 2014-8-3

编辑：Image Analyst 2014-8-3

在 MATLAB Online 中打开

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 2014-8-3

在 MATLAB Online 中打开

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!

请先登录，再进行评论。