Calculating the slope of a spectral slice (PSD)

Ok, so I have managed to create some Thompson multitaper power spectral density estimates (PSD) and with some help, I was able to also plot the mean spectrum.

What I'm trying to do now is figure out the mean of the mean spectrum , that is the ' center of gravity (COG) ,' but I want to exclude the first 1000 Hz because the power from voicing skews the COG measures quite a bit.

Based on the mean, I want to calculate, two slopes : M1 from the 1000 Hz to the mean and M2, from the mean to the end of the frequency range . I am completely lost on how to do this. The script I have so far, is as follows:

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;

h = [0.5, 0.5, 0.5];
figure1 = figure;
axes1 = axes('Parent',figure1,'YGrid','on',...
    'XTickLabel',{'0','5','10','15','20'},...
    'XTick',[0 5000 1e+004 1.5e+004 2e+004],...
    'XGrid','on');

hold(axes1, 'all');

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

grid on
title('Thompson Multitaper Power Spectral Density Estimate')
xlabel('Frequency(kHz)')
ylabel('Power/Frequency(dB/Hz)')
axis([0 22000 -150 -40])

hold off

I also want to calculate where the max peak is on the x-axis . I managed to calculate the y-value using the following:

z = max(max(meanpxx))

Which gives me z = 2.6765e-006

I then did the following to convert it into dB/Hz

z1 = 10*log10(z)

Which gives me z1 = -55.7243

How do I find the corresponding x-coordinate (aka the Frequency at the peak).

I attached a .zip with all the wav files I'm using and I attached an example plot of what I have and an example plot of what I want to create. Any help would be appreciated!

Answer 1

[Star Strider]

The max function will give you the index of the maximum value if you ask it:

[C,I] = max(A);

finds the indices of the maximum values of A, and returns them in output vector I. If there are several identical maximum values, the index of the first one found is returned.

For the slope, I would use polyfit with a first-degree polynomial.

[Phil]

Ok, the first function worked well! It gave me an I function which equals 97. This is where I'm getting stuck, it's the 97th cell, but I have a separate value for x, which is not a vector, it's just a number 44100. How do I figure out where 97 will fall on my x-axis?

[Star Strider]

Do all fa ...fi vectors have the same values and lengths? If so, then fm = fa(97) (where fm is the frequency at which the maximum value of y occurs). If they’re not all the same, I note you chose fa to plot meanpxx, so I would choose it (as I did here).

To plot it, for instance with a marker, I would do:

plot(fm, 10*log10(44100), 'pg')

to plot it with a green star. This assumes that I understand you correctly.

[Phil]

Ok! That totally worked. Thanks. I will play around with the polyfit function for a few days and see if I can make it all work. Thanks.

And yes, all the fa ... fi vectors are the same length and values.

[Star Strider]

My pleasure!

The polyfit function only accepts vector inputs, so if you want to fit 10*log10(meanpxx), I would do it as:

fit_idx = 97:length(fa);  % Index range from max(y) to the end of ‘fa’

b = polyfit(fa(fit_idx), 10*log10(meanpxx(fit_idx)), 1);
fit_meanpxx_dB = polyval(b, fa(fit_idx));

then to plot it:

plot(fa(fit_idx), fit_meanpxx_dB)

should work.

[Phil]

Hey, that does work! That's really helpful. I appreciate the help.

I plotted a second line back to the start using the following:

fit_idx2 = length(fa:0):I; %I used I here because it represents the mean and I can run multiple files without having to reenter the mean

c = polyfit(fa(fit_idx2), 10*log10(meanpxx(fit_idx2)), 1);
fit_meanpxx_dB2 = polyval(c, fa(fit_idx2));

plot(fa(fit_idx2), fit_meanpxx_dB2, 'k-', 'LineWidth', 2);

The two lines don't touch perfectly at the mean, is that normal or expected?

I have one final question, is there a way to calculate the slope in a raw number or to have MATLAB tell me what it is? For example: slope = -3?

[Phil]

Ok, I think I found it out! I used:

x = polyfit(fa(fit_idx), fit_meanpxx_dB,1)

Thanks for all the help! Seriously. Amazing.

[Star Strider]

My pleasure!

I apologise for the delay in replying.

The slope will be c(1) as calculated by your polyfit call. (The intercept will be c(2).)

I don’t quite understand ‘the two lines don’t touch perfectly at the mean’. A linear least-squares line will pass through the arithmetic mean of the x and the arithmetic mean of the y data. If you want to plot the data from the y-axis to the data at index 97, you can simply connect the y-intercept to the data at index 97, or if you want to calculate the plot:

fit_2 = polyval(c, fa(1:97));
plot(fa(1:97), fit_2)

I’m using fit_2 here because I’m not sure what you called your second vector. That should produce a single line with no gaps or other problems. Alternatively, you can plot the line over the entire range of fa with one call with this:

fit_meanpxx_dB = polyval(c, fa)
plot(fa, fit_meanpxx_dB2, 'k-', 'LineWidth', 2);

which is what I would do.

[Phil]

Ok, here is the whole script I have for plotting the slopes:

[C, I] = max(meanpxx);
C1 = 10*log10(C);
fm = fa(I);

fit_idx = I:length(fa);  % Index range from max(y) to the end of ‘fa’

b = polyfit(fa(fit_idx), 10*log10(meanpxx(fit_idx)), 1);
fit_meanpxx_dB = polyval(b, fa(fit_idx));

plot(fa(fit_idx), fit_meanpxx_dB, 'k-', 'LineWidth', 2);

fit_idx2 = length(fa:0):I;  % Index range from max(y) to the end of ‘fa’

c = polyfit(fa(fit_idx2), 10*log10(meanpxx(fit_idx2)), 1);
fit_meanpxx_dB2 = polyval(c, fa(fit_idx2));

plot(fa(fit_idx2), fit_meanpxx_dB2, 'k-', 'LineWidth', 2);

What I mean is that two slopes do not meet at the mean, the y-axis point that they both arrive at when they get to 97 is different. Is that to be expected?

I attached an image to show you what I mean, but you have to zoom rather close to see it.

Thanks again for the tremendous help and don't worry about the late replies, you have helped me a ton!

[Star Strider]

‘What I mean is that two slopes do not meet at the mean, the y-axis point that the both arrive at when they get to 97 is different. Is that to be expected?’

Yes. It is definitely to be expected. I didn’t realise you were plotting two different lines. They will not necessarily meet at the mean, because they are regressions on two different sets of data. That’s simply the nature of regression. Regression coefficients are estimates.

My pleasure! It was actually fun for me.

[Phil]

Is there some way to plot one continous line from the start, through the mean and then to the end? Of course, I don't want a continuous slope, I want to separate slopes like I have now! If not, then ok... that was super helpful though! Still a major newbie, so learning MATLAB as my first programming language takes some time. Of course, there is tons of great help available, which also makes the process enjoyable for me!

[Star Strider]

It is possible to plot a continuous line from the start through the end, probably with polyfit and a much higher order polynomial, but it would not yield the same information, and the lines would likely not be straight. You’re likely better off plotting two linear fits.

You could force both lines through the mean easily enough (essentially fixing it as the y-intercept), but that could mean sacrificing some accuracy with the slopes. It would involve linear regression — easy enough in MATLAB — but not with polyfit. If you need confidence intervals on the slope (the only parameter you would be estimating if you forced the line through the mean), it would be easiest to do with the Statistics Toolbox (particularly the regress function) if you have access to it.

[Phil]

Ok, the paper I was reading and trying to learn to do this from said they used a regression line. If that is the case, then I should certainly use the regress function! Would it be essentially the same input except I change the polyfit to regress?

Also, I am looking this up, but if you happen to answer before I figure it out, I want it to only go from 1000 Hz on the x-axis to 11025 Hz. Thanks for all the help with this!

[Star Strider]

Essentially, yes. It is just a bit more complicated than polyfit. I’d take time to explain it tonight, but I just deleted over 150 bot-posted spam messages from ‘jib ji’ over the last hour or so, and it’s posting faster than I can delete them. I’m tired (GMT-6 here).

That said, it’s probably relatively easy to set up regress to do what you want. After that, you can probably set up a function that would take your variables as arguments, call regress, and produce the values you want.

[Phil]

Sounds interesting, I will tinker with it for a bit before I go to bed.

I dunno why bots spam on here, one posted on a thread of mine before too.