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

Question

Phil 2014-8-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/146231-how-can-i-use-regress-function-to-fit-a-line-to-an-irregular-line

评论： Star Strider 2014-8-5

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2014-8-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/146231-how-can-i-use-regress-function-to-fit-a-line-to-an-irregular-line#answer_147372

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Phil 2014-8-4

编辑：Phil 2014-8-4

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

在 MATLAB Online 中打开

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.