WOW long gap to answer this question, maybe if I have a Einstein brain I can understand the image, what is not my case lol.
I haved touched in this algorithm in a long time ago, and I remember that was a surprised for me that almost no one around the Web talks about the Algorithm from Morita Naotaka, I tryed seach his original paper but I cant found online, I think that it is in Japanese, well I dont know ...
So I learned about it looking at among number of patents that make references to the Morita Naotaka, and I freeze when I see a lot of patents around 2000 years that use the same code.
As far as I know this algorithm was derived from TDHS, has advantages about TDHS, this does the job using just Period size Buffer. And believe me it's as easy to code as TDHS.
So lets go, take a look at the Yamaha Patent number US8085953 this is basicly the PICOLA from Morita Naotaka...
Here one image from this Patent that is infinity better to understand how this works to compress :
Looks like one TDHS using a triangular window, but take a look how it are cutting wave A and wave B in the exact Period size(Lp), then they are crossfaded A + B. this image seems to be an compression using rate = 0.6 factor(alpha), L can be determined using:
$$ \begin{align}
L & = \operatorname{round}\left(\frac{L_p}{\alpha-1}\right)
\end{align} $$
And here how its works to expansion:
Now L can be determined using:
$$ \begin{align}
L & = \operatorname{round}\left(L_p\frac{\alpha}{1-\alpha}\right)
\end{align} $$
Time ago I show a basic code to compress based in TDHS here, now lets try a litle bit more complex and complete example to expand (it will slowdown), so based in this images/equations, here my PICOLA Algorithm:
alpha=0.5;
f=735;
Fs=44100;
signal= 0.9*sin(2*pi*f/Fs*(1:44100)); signal = signal';
period= 60; % is it (Fs/F)
[signal,Fs]=audioread('ederwander.wav');
nsamples = length(signal);
out=zeros(1,floor(nsamples/alpha));
inptr=1;
outptr=1;
minF = 50;
maxF = 900;
minP = floor(Fs/maxF);
maxP = floor(Fs/minF);
winsize = 1024;
while ( inptr+winsize <= nsamples )
%chunk used to find the period, its need be >= MaxP
windowedSignal=signal(inptr:inptr+winsize);
%Basic Autocorrelation to find the period
Maximum=-Inf(1);
for P =minP:maxP
ac = sum(windowedSignal(1:(winsize)-P) .* windowedSignal(P:(winsize -1)));
if ac > Maximum
period=P;
Maximum=ac;
end
end
period = period-1;
%split all how shown in the patent picture
waveA = signal(inptr:inptr+period-1);
waveB = signal(inptr+period:(inptr+period*2-1));
CrossfadeAB = (waveA .* linspace(0,1,period)') + (waveB .* (1-linspace(0,1,period))');
L = round(period*alpha/(1-alpha));
waveB_T0 = signal(inptr+period:inptr+L);
out(outptr:outptr+period-1)=waveA;
out(outptr+period:outptr+period*2-1)=CrossfadeAB';
out(outptr+period*2:outptr+period+L)=waveB_T0';
%increment input and output pointers
inptr = inptr + L;
outptr = outptr + L + period;
end
%not best choice, just append the end to match the out size, this will click, maybe need crossfade here
if outptr < nsamples/alpha
fim=floor(nsamples/alpha)-outptr;
out(outptr:floor(nsamples/alpha))=signal(nsamples-fim:nsamples);
end
plot(out);
sound(out,Fs);
I think that standard algorithm of PICOLA just get nice results to 0.5<=alpha<=2
And one more time it's can be considered a proof of concept to show how its can be done!
