MATLAB: Help solving for and/or while loop.

discretizationfor loopif conditionwhile loop

Hi all,
I have a loop which gives me the frequency increments as a final result. These increments are obtained by dividing are under the function into equal area segments. I will copy paste the loop here to make it more clear:
Hs=3;
Tz=8;
Tp=1.408*Tz;
wp=2*pi/Tp;
u=20; % initial number of harmonics w_axis=3;
spacing_w=5000;
w=linspace(0,w_axis,spacing_w);
Area_under_spectrum =@(w) (5/16*Hs^2.*(wp^4./w.^5).*exp(-5/4.*(w./wp).^(-4)));
Area_total=integral(Area_under_spectrum, 0, w_axis);
area=Area_total/u;
for index=2:u+1;
    area_under_frequency_bandwidth=@(w) integral(Area_under_spectrum,0,w);
    g_function=@(w) area_under_frequency_bandwidth(w)-(index-1)*area;
    w(index)=fzero(g_function,[1e-10,w_axis]);
    frequency(index)=w(index);
    frequency_values=frequency(1,1:index);
    frequency_values(1)=0;
    delta_omega(index)=abs(frequency_values(index)-frequency_values(index-1));
end
Now, the problem is that at very low or very high frequencies the frequency increments (delta_omega) is relatively high and I want to limit it to some maximum value. Therefore, at very low or very high frequencies I want the area of a single segment to be obtained in a following manner: if a particular segment requires a frequency increment greater than the maximum I specify I want the segment area to be continuously halved until the required increment is less than the maximum.
At the end I want to have the frequency increments for all segments.
Please help me guys as I cannot solve it ;/
Thanks in advance!

Best Answer

In my view this is mostly a software design issue. You probably don't see how to do it because of the way your code is organized in one big thing.
One thing that isn't clear to me is whether you intend to divide the interval [w(i-1),w(i)] evenly to meet the maximum criterion or divide it into a power of 2 segments with equal area. Either can be done easily enough, but the code is different, obviously.
I would approach this with two basic ideas:
  1. Take the principles of the code above and write function that computes the next w value with a given w value and the target area for a segment. You probably want this function to be something like wnew = equalAreas(fun,wold,area); In that case the integral call would be based on integral(f,wold,wnew) rather than integral(f,0,wnew), and the function that goes into fzero would subtract the area variable rather than (i - 1)*area.
  2. Subdivide the resulting interval if necessary. If the interval is to be evenly divided, you'll need pow2(ceil(log2((w(i+1)-w(i))/maximum))) subintervals (just needs a nested loop). If, rather, the area is to be subdivided evenly, you'll use the equalAreas function with the third input being area/2^m for m = 1, 2, 3, ... until the maximum delta is less than the max.
Related Question