How can I use Regress Function to fit a line to an irregular line?

I've been advised to use the regress function to plot a fit for the mean spectral slice (PSD) I have created. Basically, I want to fit a line to an irregular line shape . However, I would like it to pivot at the mean , so basically, there are two slopes , one going upwards towards the mean, and then one going downwards away from the mean.

With that being said, I tried the following:

fitzone = 30:333; %indexes the vectors we want from the matrix

k = regress(fa(fitzone), 10*log10(meanpxx(fitzone)));

fitzone is there because I only want to use a portion of the vector .

Sufficed to say, I'm lost. The explanation and example with carsmall is not entirely clear to me! I'm not sure what I should be putting in for my regress(y, x). The primary data I want to fit the line to is meanpxx , but since it has gone through a multitaper power spectral density estimate I have to multiple it my 10*log10 in order to get the actual spectral shape. Below, I have pasted my entire code before the regress function and attached some sound files and an example plot (without the regression line). Any help is appreciated.

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([1000 11025 -150 -40])

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

Thanks for any help!

Answer 1

[Star Strider]

This is a preliminary version of my code. I use created data but data that are similar to yours. Note that since I’m forcing both lines to meet at the mean, the slopes are not as they would be if we estimated the intercepts as well. I can probably make a function out of it that accepts your fa and 10*log10(meanpxx) as arguments, and returns the slopes and the calculated lines as output. I’ll wait until you have a chance to run this and see if it meets your requirements.

The logic is straightforward. I put the x and y coordinates of the mean at (0,0), divided the data into two segments to the ‘left’ and ‘right’ of the mean, then computed the resulting slopes, forcing the y-intercept to be zero in both regressions. I then re-centred the calculated regression lines to overplot your data:

x = linspace(0,25);
y = [3.*x(1:20)-80  -2.*x(21:end)-55] + randn(1,100);

Lx = length(x);
mnix = find(y == max(y));               % Index of mean here
xmn = x(mnix);                          % X-Value of mean
ymn = y(mnix);                          % Y-Value of mean

xmc = x-xmn;                            % Shift X to put ‘mean’ at x=0
ymc = y-ymn;                            % Shift Y to put ‘mean’ at y=0

ix_rng1 = 1:mnix;                       % First segment
Lx1 = length(ix_rng1);
M1 = [xmc(ix_rng1)']\[ymc(ix_rng1)']    % Slope of First Segment

ix_rng2 = mnix:Lx-mnix;
Lx2 = length(ix_rng2);
M2 = [xmc(ix_rng2)']\[ymc(ix_rng2)']    % Slope of First Segment

RL1 = [x(ix_rng1)']*M1 + ymn-M1*x(mnix);
RL2 = [x(mnix:Lx)']*M2 + ymn-M2*x(mnix);

figure(1)
plot(x, y)
hold on
plot(x(ix_rng1), RL1, 'r') 
plot(x(mnix:Lx), RL2, 'r') 
hold off
grid

I apologise for the delay, but I’m still tired after last night’s spambot battle.

[Phil]

This is really awesome and don't worry at all about late replies. You are helping me out of your own kindness, so there is no need to feel a rush on my behalf.

This is almost identical to what I want, yes... what I was thinking of doing is slightly different and I have tried to no avail. What I was looking to do is let the long slope (the one with a negative slope) to just plot as it normally would and then force the second (short positive slope) to end at the start point of the other.

This shouldn't lose as much information and really only forces the slope of the shorter (left most) line. Does that make sense?

Anyway, I'm trying to change your code to do that right now, but I'm not sure I can figure it out!

One place I keep getting confused is the y input. Since I don't have two vectors as such. I have meanpxx which is a vector and fs, which is the frequency, i.e. a single cell, not a vector, so I always get confused on what to put in y coordinates.

[Star Strider]

That would not necessarily result in the two lines meeting at the mean, but if that’s not a problem, then this may do what you want:

x = linspace(0,25);
y = [3.*x(1:20)-80  -2.*x(21:end)-55] + randn(1,100);

Lx = length(x);
mnix = find(y == max(y));               % Index of mean here
xmn = x(mnix);                          % X-Value of mean
ymn = y(mnix);                          % Y-Value of mean

xmc = x-xmn;                            % Shift X to put ‘mean’ at x=0
ymc = y-ymn;                            % Shift Y to put ‘mean’ at y=0

ix_rng2 = mnix:Lx;
Lx2 = length(ix_rng2);
B2 = [ones(Lx2,1) x(ix_rng2)']\[y(ix_rng2)']    % Slope and Intercept of Second Segment
RL2 = [ones(Lx2,1)  x(mnix:Lx)']*B2;    % Fit Second Segment

ix_rng1 = 1:mnix;                       % First segment
Lx1 = length(ix_rng1);
M1 = [xmc(ix_rng1)']\[y(ix_rng1)'-RL2(1)]    % Slope of First Segment

RL1 = [x(ix_rng1)']*M1 + ymn-M1*x(mnix);
RL2 = [ones(Lx2,1)  x(mnix:Lx)']*B2;

figure(1)
plot(x, y)
hold on
plot(x(ix_rng1), RL1, 'r') 
plot(x(mnix:Lx), RL2, 'r') 
hold off
grid

It computes the slope and intercept for the negative-slope data, then subtracts the value at the mean for that regression at the x-value of the mean, and fits the positive-slope segment, forcing the intercept of the positive-slope value to be zero, just as before but with the mean at (0,0) being defined as the negative-slope regression value at that x-value.

[Phil]

This is a bit different. I think my major confusion is figuring out how to change it for my data. However, what I meant to say is I wanted to see if I can force the little slope (M1) to end exactly where M2 begins without forcing M2 to be anything except starting at the x-coordinates of the mean .

[Star Strider]

That’s difficult because of the nature of the regression and the way they’re plotted. The will take on the same values at the x-coordinate of the mean because I forced the positive-slope line to be the same as the value of the negative-slope line at that point. They won’t plot precisely (they’ll either not meet or they’ll cross) because that’s the nature of the way they’re defined and the way plot works.

With:

ix_rng2 = mnix:Lx;                      % Second Segment

and:

ix_rng1 = 1:mnix-2;                     % First Segment

there will be a slight gap between them at the mean, otherwise they would cross. Experiment with ix_rng1 and ix_rng2 to get the result you want.

I can put this in a function form to make it essentially a one-liner for you if you’re happy with the results it currently gives.

[Phil]

Thanks! You've helped me out a ton. I can't thank you enough.

[Star Strider]

My pleasure!

As scientists, we’re all in this together!

Answer 2

[dpb]

Look at John d'Errico's package

http://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling

[Phil]

This is nice and I would use it, but unfortunately, I don't think my license for MATLAB has all of these functions!