Smoothing piecewise sin functions

I wrote a function that pieces together two different sine functions. The positive curves are always sin(x) but the negative curves are sin(mult.*x).

I tried spaps/csaps/etc and spaps with tol=1 seems to give me the smoothest result, except the amplitude of the curve decreases from the desired y=-1:1. I tried adjusting the tol to different values but that didn't quite work. Is there some other function/parameter I could use to preserve both the smoothness at tol=1 and the amplitude?

The code is below:

    function [t, values]=graphic(mult);
  clf

x=0:pi./32:4.*pi;
y=[sin(x).' x.'];
%y=y(y>=0);
for m=1:length(x)
    if y(m) >=0
        y(m)=y(m);
    else
        y(m)=0;
    end
end

xx=0:pi./32:4.*pi;
yy=[sin(mult.*xx).' xx.'];
%yy=yy(yy<0);
for n=1:length(xx)
    if yy(n) < 0
        yy(n)=yy(n);
    else
        yy(n)=0;
    end
end

Y=([yy y]);
SET1=find(Y(:,3)~=0);
SET2=find(Y(:,1)~=0);
tol=1e-15;
ZERO1=find(Y(SET1,3)<tol);
ZERO2=find(Y(SET2(1):end,1)==0,1);
Inter=[Y(SET1(1:ZERO1(1)),3); Y(SET2(1:ZERO2(1)),1)];
INTER=[0; Inter; Inter];

t=1:length(INTER);
%plot(t,INTER);
%y2=smooth(INTER,'sgolay')
% p=polyfit(t,INTER.',6);
% y2=polyval(p,t);
%plot(t,y2);
[sp,values]=spaps(t,INTER,1);
fnplt(sp);
set(gca,'XTick',0:25:length(INTER));
end

[dpb]

Don't think your above code works or if does run doesn't do what you think/want...

"The Matlab way"... :)

x=0:pi/32:4*pi;
y=sin(x)';
y(y<0)=0;     % above you've already built a 2-col matrix but
y=[y x'];     % address it only w/ one subscript.

Not sure what the purpose of the x in the y array is for???

Ditto comments above for next section plus your 'xx' is identical to 'x'; why not just use the copy you already have?

yy=sin(mult.*x)';  % do the fixup before the concatenation
yy(yy>=0)=0;
yy=[yy x'];          % again, why the concatenation into yy

why don't you just write

x=[0:pi/32:4*pi]';         % Start w/ a column if that's what want...
y=sin(x);                  % then y is already column, too.
iy=y<0;                    % logical vector of the y<0 locations
y(iy)=0;                   % clear them
y(~iy)=sin(mult*x(~iy));   % fill in the alternate locations

What's mult? Not clear how to answer the rest of the question w/o that minor detail...

[meihua]

For future reference, I used the spline fitting toolbox function csapi to approximate the same kind of line created by the code above.

Answer 1

[Image Analyst]

Your code doesn't run, and you didn't upload a diagram, so it's hard to imagine what it looks like. But I imagine stretches of sine waves with different periods (frequencies) and where they meet you are going to have a huge step in intensity.

To smooth the steps out, you can use Fourier filtering (or use conv2()), or you can use a Savitzky-Golay filter. A sine can be pretty well fit by a 3rd or 5th order polynomial in the short run, so if you have the Signal Processing Toolbox, you can use sgolay() with an order 3 or 5 and it should do a pretty good job of leaving the sine wave parts along but smoothing over the abrupt steps where the sines of different frequencies meet.

[meihua]

Edited my code to include the function line. Mult is the parameter that modifies the negative sine curve. Already tried sgolay and other smoothing functions, and it was too abrupt at the meeting points on the x axis. That's why I went for the spline method.

[dpb]

Well, that was already clear--the question is, what's typical value in comparison to the alternate case where

mult=1

???

How disparate the sections are depends upon that mismatch magnitude.

[meihua]

Mult=2,4, or 8.

[Image Analyst]

The function line is wrong but after I fixed that it complained about spaps(). You didn't list any other toolboxes in the Product box above, so it's probably some function you custom wrote that you didn't supply. Bottom line, I still can't run it.

[meihua]

Yes...it's a custom function...the function line written exactly is

Function [t,values]=graphic(mult)
 Insert above script
End

Then type graphic(2) in the command window...

[Image Analyst]

Huh? You give me more non-MATLAB code? Why??? That's not even legitimate MATLAB code. What about spaps()? That's what I had a problem with, not graphic().

[dpb]

>> which spaps
C:\ML_R2012b\toolbox\curvefit\splines\spaps.m

>> help spaps
spaps Smoothing spline.

    [SP,VALUES] = spaps(X,Y,TOL)  returns the B-form and, if asked, the
    values at X, of a cubic smoothing spline  f  to the data X, Y.

Despite that, I think OP will be better off to generate a continuous function instead by varying his frequency instead.

Overlooking his problems w/ posting working code snippets, just generate sin(x) from [0-pi] then -sin(2x) from [pi-2pi] and you have his starting point for the case mult=2 which has the discontinuity in omega that he's trying to smoosh out after the fact. Better to not introduce it to start with. (imo, $0.02, ymmv, etc., etc., etc., ... :) )

[Image Analyst]

I added it to the Product list for him. I don't have that toolbox.

[meihua]

Look, I edited my original post to include the entirety of my script as is. spaps came in my installed version of Matlab so I wasn't aware you didn't have the toolbox.

Answer 2

[dpb]

Instead of trying to mung on the result of mismatching frequencies caused by the discontinuous mult value as you have, redefine it to be a continuous function w/ a rapid change.

A sigmoid or the like would work to remove the discontinuity but still allow a rapid change from one to the other frequency and vice versa.

ADDENDUM: Example for first cycle of switch from wx=1x to wx=2x as a function of changing the intensity of the sigmoid. Try the following from command line...

x=[0:pi/128:2*pi]';     % need better resolution than your original...
N=[0.1 0.5 1 2 5 10]';  % a set of weights for the sigmoid--read code
M=2;                    % the target multiplier value
s=zeros(length(x),length(N));  % preallocate so can plot all together later
for i=1:length(N)
  n=N(i);
  w=(M-1)*(tanh(n*(x-pi))+1)/2+1;  % the omega weighted vector
  s(:,i)=sin(w.*x);                % and the sin resulting
end
plot(x,s)
legend(num2str(N))

You'll observe that as you make the transition sharper and sharper the resulting discontinuity approaches that of the limiting case--I had envisioned you could get by "more better" than it actually turns out, but certainly the way to generate a smooth transition is something on this order rather than trying to fix up the abrupt transition as currently trying to do.

Doing the obverse transition back at 2pi is left as "exercise for the student" -- enjoy :)

OBTW, as did in previous comment, it's worthwhile to plot the weight function and/or just the sigmoid as function of the various parameters to compare the effects w/ the result.

ADDENDUM 2:

On the comment of modifying the pi:2pi regions only, change the offset bias from the above pi to something >pi: that will be the breakpoint. Plotting the sigmoid for various choices of parameters is best way to see what's going on and what affects what...

[meihua]

Can you give me an example of this function? Ideally, I want to progressively decrease the width of the negative curves only.

[Image Analyst]

Can you first give us an example of code that actually works so we can generate your data?

[dpb]

Sample sigmoid scaled to range 0-1.

x=[0:pi/32:4*pi]';
xx=[-fliplr(x) x];
plot(xx,(tanh(xx)+1)/2)

To use, since sin(x)<0 for the range of multiples of [pi < x < 2pi], shift the origin to pi and scale the argument to create as narrow a band as desired. Then fliplr at opposite end.

Build the waveform for mult over the interval 0-2pi and use the .* operators to evaluate simultaneously over that range. repmat that result as often as needed.

[meihua]

I do need to repeat the curve from 0:2pi but now what I'm seeing is a discontinuity at 2pi instead of pi.

[dpb]

See sample code in Answer section -- you'll have to build the sigmoid at the two points--pi and 2pi. I demonstrated it at pi; as noted you'll want the reverse sense at the other end to go from w=Mx back to w=1x. I left this specific exercise "for the student" -- show what you've done to date. Did you look at the plotted sample case and follow how it works as the pattern?

[meihua]

x=[0:pi/32:4*pi]';   
m=.5*(tanh(.5*(pi+x))+.5*x-tanh(.5*(2*pi-x)))/2+1
plot(x,sin(m.*x));

I don't think this is what you meant. I tried fiddling around with it a lot, and I kind of understand what you're saying to do--I just don't know what I should be manipulating to get that "obverse transition back at 2pi."

[dpb]

Couple of ways...

a)

w2=(M-1)*(tanh(n*(x(floor(end/2):end)-3*pi))+1)/2+2;

where now

x=0:4*pi

interval so w2 is covering 2pi-4pi and w is as before except for 0-2pi or x(1:round(end/2))

b) x as above over the extended range, then create w as

w=[w;flipud(w)];

May need to fixup length(w) by an "off-by-one" error if length(x) is odd the length of doubled w will be one too great so use (end-1:-1:1) for the added portion to match the length and not duplicate the middle point.

[meihua]

n=3;
M=2;
%x=[0:pi/32:4*pi]';
x=linspace(0,4*pi,150)';
xx=[-fliplr(x) x];
w=(M-1)*(tanh(n*(x(1:round(end/2))-pi))+1)/2+1;
w2=(M-1)*(tanh(n*(x(round(end/2)+1:end)-3*pi))+1)/2+2;
%m=.5*(tanh(.5*(pi+x))+.5*x-tanh(.5*(x+2*pi)))/2+1;
W=cat(1,w,w2);
plot(x,sin(W.*x));
xlim([0 4*pi]);

This is what I got by plugging in w2. But the graph looks very strange, and manipulating the two equations doesn't produce much effect...