function MKS_04a
% MKS_04a - FÁZOVÝ VOKODÉR. Výuková ukázka v MATLABu bez dalších toolboxů.
% Bez korekce očekávaného fázového posunu a bez převodu rozdílů do hlavního intervalu fáze.
% Nepřekrývající se syntetické rámce mohou způsobit mezery a slyšitelné artefakty.
% Soubory why11.wav, Tom.wav a trumpet.wav uložte vedle skriptu (nebo použijte LOAD WAV).
root=fileparts(mfilename('fullpath'));
S.fs=22050; S.x=[]; S.y=[]; S.player=[]; S.name='';
fig=figure('Name','MKS 04a - PHASE VOCODER','NumberTitle','off', ...
    'MenuBar','none','ToolBar','none','Resize','off', ...
    'Color',[.95 .96 .97],'Position',[110 65 1100 770], ...
    'CloseRequestFcn',@closeFigure);
bg=get(fig,'Color');
uicontrol(fig,'Style','text','String','MKS 04a - PHASE VOCODER', ...
    'Position',[230 725 640 32],'FontSize',19,'FontWeight','bold','BackgroundColor',bg);
uicontrol(fig,'Style','text','String','SOURCE','Position',[35 683 90 22], ...
    'FontWeight','bold','FontSize',11,'BackgroundColor',bg);
source=uicontrol(fig,'Style','popupmenu', ...
    'String',{'why11.wav','Tom.wav','trumpet.wav','SINE 440 Hz','LOAD WAV ...'}, ...
    'Position',[128 682 225 28],'FontSize',11,'Callback',@sourceChanged);
fileText=uicontrol(fig,'Style','text','String','', ...
    'Position',[370 683 685 23],'HorizontalAlignment','left','BackgroundColor',bg);

panel=uipanel(fig,'Title','ANALYSIS / SYNTHESIS','FontWeight','bold', ...
    'Units','pixels','Position',[30 552 1040 117],'BackgroundColor',bg);
makeLabel(panel,'WINDOW LENGTH N',30,65,170);
winPop=uicontrol(panel,'Style','popupmenu','String',{'32','64','128','256','512','1024','2048','4096'}, ...
    'Value',6,'Position',[30 37 170 27],'Callback',@parametersChanged);
winMsText=makeLabel(panel,'',30,12,185);
makeLabel(panel,'OVERLAP',230,65,160);
overlapPop=uicontrol(panel,'Style','popupmenu','String',{'0 %','25 %','50 %','75 %','87.5 %'}, ...
    'Value',4,'Position',[230 37 160 27],'Callback',@parametersChanged);
makeLabel(panel,'TIME FACTOR R',430,65,230);
ratioSl=uicontrol(panel,'Style','slider','Min',.2,'Max',5,'Value',1, ...
    'Position',[430 42 215 20],'Callback',@parametersChanged);
ratioText=makeLabel(panel,'R = 1.00',650,39,110);
makeLabel(panel,'PITCH SHIFT (semitones)',775,65,225);
pitchPop=uicontrol(panel,'Style','popupmenu', ...
    'String',arrayfun(@(v)sprintf('%+d',v),-12:12,'UniformOutput',false), ...
    'Value',13,'Position',[775 37 210 27],'Callback',@parametersChanged);
info=uicontrol(fig,'Style','text','String','', ...
    'Position',[45 515 1010 25],'HorizontalAlignment','center', ...
    'FontSize',10,'BackgroundColor',bg);
axAnalysis=axes('Parent',fig,'Units','pixels','Position',[70 425 960 66]);
axSynthesis=axes('Parent',fig,'Units','pixels','Position',[70 329 960 66]);
axInput=axes('Parent',fig,'Units','pixels','Position',[80 116 430 170]);
axOutput=axes('Parent',fig,'Units','pixels','Position',[605 116 430 170]);
uicontrol(fig,'Style','pushbutton','String','PLAY ORIGINAL', ...
    'Position',[48 53 160 42],'FontSize',11,'Callback',@playOriginal);
uicontrol(fig,'Style','pushbutton','String','PROCESS', ...
    'Position',[225 53 160 42],'FontSize',11,'FontWeight','bold','Callback',@processAudio);
uicontrol(fig,'Style','pushbutton','String','PLAY RESULT', ...
    'Position',[402 53 160 42],'FontSize',11,'Callback',@playResult);
uicontrol(fig,'Style','pushbutton','String','STOP', ...
    'Position',[579 53 125 42],'FontSize',11,'Callback',@stopAudio);
uicontrol(fig,'Style','pushbutton','String','SAVE WAV', ...
    'Position',[721 53 125 42],'FontSize',11,'Callback',@saveWav);
status=uicontrol(fig,'Style','text','String','', ...
    'Position',[35 13 1030 28],'HorizontalAlignment','left', ...
    'FontSize',10,'BackgroundColor',bg);
sourceChanged([],[]);

    function h=makeLabel(parent,str,x,y,w)
        h=uicontrol(parent,'Style','text','String',str,'Position',[x y w 21], ...
            'HorizontalAlignment','left','BackgroundColor',bg,'FontSize',10);
    end
    function sourceChanged(~,~)
        items=get(source,'String'); idx=get(source,'Value'); name=items{idx};
        if idx==5
            [fn,fp]=uigetfile('*.wav','Select WAV');
            if isequal(fn,0),return;end
            path=fullfile(fp,fn); name=fn;
        elseif idx==4
            S.fs=22050; t=(0:2*S.fs-1)'/S.fs;
            S.x=.65*sin(2*pi*440*t); S.name='SINE 440 Hz';
            set(fileText,'String',sprintf('Synthetic sinusoid, %.2f s',numel(S.x)/S.fs));
            invalidate(); return;
        else
            path=fullfile(root,name);
            if ~exist(path,'file')
                [fn,fp]=uigetfile('*.wav',['Find ' name]);
                if isequal(fn,0),set(status,'String',['Missing file: ' name]);return;end
                path=fullfile(fp,fn); name=fn;
            end
        end
        try
            [x,fs]=audioread(path);
            x=mean(x,2); x=x-mean(x);
            S.x=x(:); S.fs=fs; S.name=name;
            set(fileText,'String',sprintf('%s   |   %d Hz   |   %.2f s',name,fs,numel(x)/fs));
            invalidate();
        catch ME
            set(status,'String',['Cannot read audio: ' ME.message]);
        end
    end
    function parametersChanged(~,~)
        invalidate();
    end
    function invalidate()
        stopAudio(); S.y=[];
        N=2^(get(winPop,'Value')+4);
        ov=[0 .25 .5 .75 .875]; Ha=round(N*(1-ov(get(overlapPop,'Value'))));
        R=get(ratioSl,'Value'); semitones=get(pitchPop,'Value')-13;
        ratioInternal=R*2^(semitones/12); Hs=Ha;
        set(ratioText,'String',sprintf('R = %.2f',R));
        set(winMsText,'String',sprintf('%.1f ms at %.0f Hz',1000*N/S.fs,S.fs));
        set(info,'String',sprintf('N = %d    Ha = Hs = %d    frame-count ratio = %.3f',N,Ha,ratioInternal));
        drawWindows(N,Ha,ratioInternal);
        if ~isempty(S.x)
            showSpec(axInput,S.x,S.fs,'ORIGINAL');
        end
        cla(axOutput); title(axOutput,'RESULT (press PROCESS)');
        set(status,'String','Set parameters and press PROCESS. Pitch shift uses resampling after time stretch.');
    end
    function drawWindows(N,Ha,ratio)
        % Stejný krok analýzy a syntézy; mění se počet syntetických rámců.
        % Krátká ukázka obou posloupností na STEJNÉ časové ose.
        w=.5-.5*cos(2*pi*(0:N-1)'/(N-1));
        colors=lines(4); nA=8;
        nS=max(1,round(nA*ratio));
        xMax=(max(nA,nS)*Ha+N)/N;
        axs=[axAnalysis axSynthesis]; counts=[nA nS];
        captions={'ANALYSIS WINDOWS','SYNTHESIS WINDOWS'};
        for p=1:2
            ax=axs(p); cla(ax); hold(ax,'on');
            for k=0:counts(p)-1
                tt=((0:N-1)+k*Ha)/N;
                plot(ax,tt,w,'Color',colors(mod(k,4)+1,:),'LineWidth',1.1);
            end
            xlim(ax,[0 xMax]); ylim(ax,[0 1.12]);
            ylabel(ax,'Amplitude');
            title(ax,sprintf('%s | hop = %d | %d frames (illustration)', ...
                captions{p},Ha,counts(p)));
            grid(ax,'on'); hold(ax,'off');
            if p==1
                set(ax,'XTickLabel',[]);
            else
                xlabel(ax,'Time (in window lengths N)');
            end
        end
    end

    function processAudio(~,~)
        if isempty(S.x),return;end
        stopAudio();
        N=2^(get(winPop,'Value')+4);
        ov=[0 .25 .5 .75 .875]; Ha=round(N*(1-ov(get(overlapPop,'Value'))));
        R=get(ratioSl,'Value'); semitones=get(pitchPop,'Value')-13;
        pitchFactor=2^(semitones/12);
        % Postup podle původního SYN04: interpolace spektrálních rámců a následné převzorkování.
        % Délka výsledku je přibližně R-násobkem délky originálu.
        internalRatio=R*pitchFactor;
        Hs=Ha; % Pevný krok syntézy; mění se počet rámců.
        % Nulový překryv analýzy zůstává dostupný pro demonstraci.
        set(status,'String','Processing ...');drawnow;
        try
            z=phaseVocoderStretch(S.x,internalRatio,N,Ha);
            target=max(1,round(numel(S.x)*R));
            if semitones~=0
                % Interpolace realizuje převzorkování bez dalších toolboxů.
                query=1+(0:target-1)'*pitchFactor;
                z=interp1((1:numel(z))',z,query,'linear',0);
            end
            S.y=fitToLength(z,target);
            showSpec(axOutput,S.y,S.fs,'RESULT');
            if Ha>=N
                warningText=' | WARNING: zero overlap; artifacts possible';
            else
                warningText='';
            end
            set(status,'String',sprintf('Input %.3f s -> output %.3f s | R = %.2f | pitch %+d st%s', ...
                numel(S.x)/S.fs,numel(S.y)/S.fs,R,semitones,warningText));
        catch ME
            S.y=[];set(status,'String',['Processing error: ' ME.message]);
        end
    end
    function playOriginal(~,~),playAudio(S.x);end
    function playResult(~,~)
        if isempty(S.y),set(status,'String','Press PROCESS first.');return;end
        playAudio(S.y);
    end
    function playAudio(x)
        if isempty(x),return;end
        stopAudio();p=max(abs(x));
        if p>0,x=.95*x/p;end % Přehrávání s normalizací podobnou soundsc; data se nemění
        try
            S.player=audioplayer(x,S.fs);play(S.player);
        catch ME
            set(status,'String',['Playback error: ' ME.message]);
        end
    end
    function stopAudio(~,~)
        if ~isempty(S.player)
            try,stop(S.player);catch,end
            S.player=[];
        end
    end
    function saveWav(~,~)
        if isempty(S.y),set(status,'String','Press PROCESS first.');return;end
        [fn,fp]=uiputfile('*.wav','Save phase-vocoder result','MKS_04a_result.wav');
        if isequal(fn,0),return;end
        y=S.y;pk=max(abs(y));if pk>1,y=.99*y/pk;end
        audiowrite(fullfile(fp,fn),y,S.fs);
        set(status,'String',['Saved: ' fullfile(fp,fn)]);
    end
    function closeFigure(~,~),stopAudio();delete(fig);end
end


function y=fitToLength(x,n)
x=x(:);y=zeros(n,1);k=min(n,numel(x));y(1:k)=x(1:k);
end

function showSpec(ax,x,fs,ttl)
x=x(:);N=1024;hop=256;
if isempty(x),cla(ax);return;end
if numel(x)<N,x(end+1:N)=0;end
K=1+floor((numel(x)-N)/hop);
% Omezíme počet zobrazených rámců, aby rozhraní zůstalo svižné.
step=max(1,ceil(K/350));frames=1:step:K;
A=zeros(N/2+1,numel(frames));w=.5-.5*cos(2*pi*(0:N-1)'/N);
for j=1:numel(frames)
    idx=(1:N)+(frames(j)-1)*hop;
    F=fft(x(idx).*w);
    A(:,j)=20*log10(abs(F(1:N/2+1))+1e-6);
end
imagesc(ax,(frames-1)*hop/fs,(0:N/2)*fs/N,A);
axis(ax,'xy');ylim(ax,[0 min(5000,fs/2)]);
xlabel(ax,'Time (s)');ylabel(ax,'Frequency (Hz)');title(ax,ttl);
end

function y=phaseVocoderStretch(x,ratio,N,Ha)
% FÁZOVÝ VOKODÉR: STFT -> interpolace amplitud -> akumulace fáze
% -> inverzní FFT a vážené překrývání a sčítání rámců (WOLA).
% Pevný krok analýzy i syntézy Ha; změna délky mění počet rámců.
% Výukové zjednodušení: přímé rozdíly fází (bez rozbalování fáze).
% 1. Krátkodobá Fourierova transformace (STFT)
x=x(:); L=numel(x); pad=N;
xp=[zeros(pad,1);x;zeros(pad+N,1)];
K=1+floor((numel(xp)-N)/Ha);
w=sqrt(max(0,.5-.5*cos(2*pi*(0:N-1)'/N)));
X=zeros(N,K);
for k=1:K
    idx=(1:N)+(k-1)*Ha;
    X(:,k)=fft(xp(idx).*w);
end
% Použijeme nezáporné frekvence a obnovíme komplexně sdruženou symetrii.
% 2. Interpolace amplitud na nové časové mřížce
H=N/2+1; A=abs(X(1:H,:)); P=angle(X(1:H,:));
q=0:1/ratio:(K-1);
if isempty(q),q=0;end
B=interp1(0:K-1,A.',q.','linear').';
% 3. Kumulativní součet fázových přírůstků
D = diff(P,1,2);                  % Rozdíly fází sousedních rámců
phase = zeros(H,numel(q));
phase(:,1) = P(:,1);                 % Počáteční fáze
if numel(q)>1
    new_grid = min(size(D,2),floor(q(2:end))+1);
    D_new = D(:,new_grid);          % Fázové přírůstky ve vybraných rámcích
    phase(:,2:end) = P(:,1) + cumsum(D_new,2);
end
% 4. Komplexní spektrum a sdružená symetrie
Y=zeros(N,numel(q));
positive=B.*exp(1i*phase);
positive(1,:)=real(positive(1,:));
positive(end,:)=real(positive(end,:));
Y(1:H,:)=positive;
Y(H+1:end,:)=conj(positive(end-1:-1:2,:));
% 5. Inverzní FFT a vážené překrývání a sčítání rámců
Kout=size(Y,2); outN=(Kout-1)*Ha+N;
acc=zeros(outN,1); normW=zeros(outN,1);
for k=1:Kout
    idx=(1:N)+(k-1)*Ha;
    frame=real(ifft(Y(:,k))).*w;
    acc(idx)=acc(idx)+frame;
    normW(idx)=normW(idx)+w.^2;
end
valid=normW>1e-10; acc(valid)=acc(valid)./normW(valid);
start=pad+1; expected=max(1,round(L*ratio));
y=zeros(expected,1);
last=min(outN,start+expected-1);
if last>=start,y(1:last-start+1)=acc(start:last);end
end