ksstringexpofn.m

非常好的数字处理教程

    end
    y = zeros(delay_length,1);

    excite_contents = get(ud.excite_menu,'String');
    switch excite_contents{get(ud.excite_menu,'Value')}
        case 'Pluck'
            y(1:excitepos) = amp/excitepos*[0:excitepos-1];
            y(excitepos+1:end) = amp/(1-excitepos)*([excitepos:delay_length-1]-excitepos)+amp;
        case 'Strike'
						% Y is velocity, not stace
            hammerl = floor(delay_length/5);
            if excitepos+hammerl > delay_length,
                excitepos = delay_length-hammerl-1;
            end
            y(excitepos:excitepos+hammerl) = amp;
        case 'Noise'
            sig = 0.5 - rand(1,delay_length-2);
            % lowpass filter
            [N, Wn] = ellipord(0.05,0.1,1,40);
            [B,A] = ellip(N,1,40,Wn,'low');
            sig = filter(B,A,sig);
            sig = amp*sig./max(sig);
            y(2:end-1) = sig;
        case 'Bow'
            y(1:excitepos) = amp/excitepos*[0:excitepos-1];
            y(excitepos+1:end) = amp/(1-excitepos)*([excitepos:delay_length-1]-excitepos)+amp;
            temp = resample(y,1,5);
            excitelen = size(temp,1);
            %bowshape = [zeros(excitepos,1); temp; zeros(size(y,1)-size(temp,1)-excitepos,1)];
            bowshape = y;
            % Should be based on bow velocity
            tablelength = 1000;
            finc = floor(tablelength*0.9);  % samples spent increasing pull
            fdec = tablelength - finc;  % samples spent flying back
            forcetable = [linspace(0,amp,finc)'; linspace(amp,0,fdec)'];
            y = zeros(delay_length,1);
            clear temp;
    end

    udl = y/2;
    ldl = y/2;
    

    % Now simulate
    y = zeros(numsamples,1);

    switch excite_contents{get(ud.excite_menu,'Value')}
        case 'Bow'
            contents = get(ud.model_menu,'String');
            switch contents{get(ud.model_menu,'Value')}
                case {'Ideal','Lossy'}
                    for i=1:numsamples
                        % Move one sample from ldl to udl, and vice versa
                        tempval = udl(end);
                        udl = [-loss*ldl(1); udl(1:end-1)];
                        ldl = [ldl(2:end); -loss*tempval];
                        for j = 1:length(pickups)
                            y(i) = y(i) + udl(pickups(j)) + ldl(pickups(j));
                        end
                        udl = udl + (bowshape/2)*forcetable(mod(i,tablelength)+1);
                        ldl = ldl + (bowshape/2)*forcetable(mod(i,tablelength)+1);
                    end
                case {'Loss(f)'}
                    for i=1:numsamples
                        % Move one sample from ldl to udl, and vice versa
                        tempval = udl(end);
                        udl = [-loss*ldl(1); udl(1:end-1)];
                        % Lowpass filter
                        newval = (tempval+udl(end))/2;
                        ldl = [ldl(2:end); -loss*newval];
                        for j = 1:length(pickups)
                            y(i) = y(i) + udl(pickups(j)) + ldl(pickups(j));
                        end
                        udl = udl + (bowshape/2)*forcetable(mod(i,tablelength)+1);
                        ldl = ldl + (bowshape/2)*forcetable(mod(i,tablelength)+1);
                    end
						end
				case 'Strike'
					contents = get(ud.model_menu,'String');
					switch contents{get(ud.model_menu,'Value')}
							case {'Ideal','Lossy'}
									for i=2:numsamples
											% Move one sample from ldl to udl, and vice versa
											tempval = udl(end);
											udl = [-loss*ldl(1); udl(1:end-1)];
%                         % Move temp val through allpass
% 
%                         if length(udl) < length(num)
%                             x =[udl; ldl(1:length(num)-length(udl))];
%                         else
%                             x = udl(end-length(num):end);
%                         end
%                             
%                         tempval = filter(num,den,x);
											ldl = [ldl(2:end); -loss*tempval(1)];
											% Leaky integrator
											for j = 1:length(pickups)
													y(i) = y(i) + udl(pickups(j)) + ldl(pickups(j)) + 0.99*y(i-1);
											end
									end
							case {'Loss(f)'}
									for i=2:numsamples
											% Move one sample from ldl to udl, and vice versa
											tempval = udl(end);
											udl = [-loss*ldl(1); udl(1:end-1)];
											% Lowpass filter
											newval = (tempval+udl(end))/2;
											ldl = [ldl(2:end); -loss*newval];
											% Leaky integrator
											for j = 1:length(pickups)
													y(i) = y(i) + udl(pickups(j)) + ldl(pickups(j)) + 0.99*y(i-1);
											end
									end
					end
        otherwise
            contents = get(ud.model_menu,'String');
            switch contents{get(ud.model_menu,'Value')}
                case {'Ideal','Lossy'}
                    for i=1:numsamples
                        % Move one sample from ldl to udl, and vice versa
                        tempval = udl(end);
                        udl = [-loss*ldl(1); udl(1:end-1)];
%                         % Move temp val through allpass
% 
%                         if length(udl) < length(num)
%                             x =[udl; ldl(1:length(num)-length(udl))];
%                         else
%                             x = udl(end-length(num):end);
%                         end
%                             
%                         tempval = filter(num,den,x);
                        ldl = [ldl(2:end); -loss*tempval(1)];
                        
                        for j = 1:length(pickups)
                            y(i) = y(i) + udl(pickups(j)) + ldl(pickups(j));
                        end
                    end
                case {'Loss(f)'}
                    for i=1:numsamples
                        % Move one sample from ldl to udl, and vice versa
                        tempval = udl(end);
                        udl = [-loss*ldl(1); udl(1:end-1)];
                        % Lowpass filter
                        newval = (tempval+udl(end))/2;
                        ldl = [ldl(2:end); -loss*newval];
                        for j = 1:length(pickups)
                            y(i) = y(i) + udl(pickups(j)) + ldl(pickups(j));
                        end
                    end
            end
    end

    % Taper the first 5 milliseconds

    taperlength = ceil(Fs*5/1000);
    if (size(y,1) < 2*taperlength)
        % Taper the first 5%
        taperlength = ceil(0.05*size(y,1));
    end

    for i=1:size(y,2)
        y(1:taperlength,i) =  y(1:taperlength,i).* ...
            sin(pi*[0:taperlength-1]'/(2*taperlength));
    end

    
    % Taper the last 50 milliseconds if sound is longer than 300 ms
    if size(y,1)/Fs < 0.3
        taperlength = ceil(Fs*50/1000);
        if (size(y,1) < 2*taperlength)
            % Taper the last 20%
            taperlength = ceil(0.2*size(y,1));
        end

        for i=1:size(y,2)
            y(end-taperlength+1:end,i) = y(end-taperlength+1:end,i).* ...
                cos(pi*[0:taperlength-1]'/(2*taperlength));
        end
    else
        for i=1:size(y,2)
            y(end-taperlength+1:end,i) = y(end-taperlength+1:end,i).* ...
                cos(pi*[0:taperlength-1]'/(3*taperlength));
        end
    end
    
    
    % Apply Instrument impulse response
    if isfield(ud,'impulse_inst')
        wetness = str2num(get(ud.wetness_inst,'String'))/100;
        impulse = ud.impulse_inst.data;

        if (size(y,1) >= size(impulse,1))
            for i=1:size(impulse,2)     % For each channel
                temp(:,i) = [y; zeros(size(impulse(:,i),1),1)];
                temp(:,i) = fftfilt(impulse(:,i),temp(:,i));
            end
        else
            impulse = [impulse; zeros(size(y,1),size(impulse,2))];
            for i=1:size(impulse,2)
                temp(:,i) = fftfilt(y,impulse(:,i));
            end
        end

        temp = temp*wetness;
        for i=1:size(impulse,2),
            temp(1:size(y,1),i) = temp(1:size(y,1),i) + y.*(1-wetness)./size(impulse,2);
        end
        y = temp;
    end

    
    % Apply Room impulse response
    if isfield(ud,'impulse_room')
        wetness = str2num(get(ud.wetness_room,'String'))/100;
        impulse = ud.impulse_room.data;
        clear temp;
        
        if (size(y,1) >= size(impulse,1))
            for i=1:size(impulse,2)     % For each channel
                temp(:,i) = [y; zeros(size(impulse(:,i),1),1)];
                temp(:,i) = fftfilt(impulse(:,i),temp(:,i));
            end
        else
            impulse = [impulse; zeros(size(y,1),size(impulse,2))];
            for i=1:size(impulse,2)
                temp(:,i) = fftfilt(y,impulse(:,i));
            end
        end
        
        temp = temp*wetness;
        for i=1:size(impulse,2),
            temp(1:size(y,1),i) = temp(1:size(y,1),i) + y.*(1-wetness)./size(impulse,2);
        end
        y = temp;
    end

%     % Taper the first 10 milliseconds if no impulse
%     if ~isfield(ud,'impulse_room') & ~isfield(ud,'impulse_inst')
%         taperlength = ceil(Fs*10/1000);
%         if (size(y,1) < 2*taperlength)
%             % Taper the first 5%
%             taperlength = ceil(0.05*size(y,1));
%         end
% 
%         for i=1:size(y,2)
%             y(1:taperlength,i) =  y(1:taperlength,i).* ...
%                 sin(pi*[0:taperlength-1]'/(2*taperlength));
%         end
%     end
    
    % Taper the last 50 milliseconds
%     taperlength = ceil(Fs*50/1000);
%     if (size(y,1) < 2*taperlength)
%         % Taper the last 20%
%         taperlength = ceil(0.2*size(y,1));
%     end
%     
%     for i=1:size(y,2)
%         y(end-taperlength+1:end,i) = y(end-taperlength+1:end,i).* ...
%             cos(pi*[0:taperlength-1]'/(2*taperlength));
%     end
    
    if max(abs(y)) > 0.9
        y = normalize(y);
    end
    
    ud.audiodata.data = y;
    ud.audiodata.Fs = Fs;

    
function [num,den] = designAllpass(delay)

    N = floor(delay);

    d = zeros(1,N);
    p = zeros(1,N+1);

    for k=1:N,
        for n=0:N,
            p(n+1) = (delay-N+n)/(delay-N+k+n);
        end
        d(k) = (-1)^k*nchoosek(N,k)*prod(p);
    end

    w = linspace(0,1,500).*pi;

    % Allpass filter giving fractional delay
    num = [fliplr(d) 1];
    den = [1 d];