% E02_3.m          08mar18  16mar23
% Impulse Response Computation - Audio manipulation examples
clear all;
return

%% impulse and step responses for FIR filters; convolution
h = [.8  -.44  .36  -.19];
N = 15; % input length: N=1 is enough for impulse resp. of an FIR system
 x = [1  zeros(1,N-1)]; % impulse;
 x = ones(1,N); % step;  N >= length(h);  NOTE: disregard last part of y (negative step)
y = conv(x,h,'same');  % default 'full', otherwise 'same'
stem(y);    % 'full' length = N+4-1

%%  'filter' cmd
help filter

%%  impulse and step responses for FIR & IIR filters;
N = 31;
  p = [.8  -.44  .36  -.19];    d = [1];
  p = ones(1,9)/9;    d = [1];
  p = 1 ./ [1:9];    d = [1];                        % lossy moving average
  p = (-1).^[0:8] ./ [1:9];    d = [1];              % hf emphasis (note step response)
  p = [1];    d = [1 -1];                            % accumulator
  p = [1];    d = [1 -.7];                           % lossy accumulator
%  p = [1];    d = [1  .7];                           % hf emphasis (note step response)
  p = [.8  -.44  .36  -.19];    d = [1  -.7]; 
%  p = [.8  -.44  .36  -.19];    d = [1   .7]; 
%  p = [.8  -.44  .36  -.44  .8];    d = [1];         % linear phase FIR

x = [1  zeros(1,N-1)];
y = filter(p,d,x);
k = 0:1:N-1;
subplot(211);  stem(k,y);  xlabel('Time index n'); ylabel('Impulse response')

x = ones(1,N);
y = filter(p,d,x);
subplot(212);  stem(k,y);  xlabel('Time index n'); ylabel('step response')

% fvtool(p,d)

%% filtering audio data   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf
% nam = 'loonytun.wav';     % 11000 Hz, 8 bit, 71924 samples
 nam = 'around_the_world-atc.wav';    % 22050 Hz, 16 bit, 4814784 samples
audioinfo(nam)
x = audioread(nam);
x = x(13e4+1:23e4);
[min(x), max(x)]

figure(1);
subplot(411); plot(x(3001:4000)); % note this is just .05s at 22kHz
subplot(412); plot(x(90e3+1:91e3));

 Fs = 22050; sound(x,Fs);
x_orig = x;


%% quantiz. noise and its frequency content
for Nb = 15:-1:1   %  no. of bits (+1 sign)   orig. 16 bit
Nb
xx = round(x*2^Nb)/2^Nb;  % linear quantizer (simplified)
figure(2);
    subplot(211); hist(x,256);
    subplot(212); hist(xx,256);
figure(3);
    subplot(311); plot(xx-x);
    subplot(312); plot(abs(fft(x)));
    subplot(313); plot(abs(fft(xx-x)));
if Nb<4, sound(xx,Fs), end
pause()
end


%%  additive impulse (long-tailed) noise - downsampling - changing Fs
x = x_orig;
figure(1);
d = rand(size(x));  d(d<.97) = 0;  d(d>=.97) = 1;  % noise prob. = 3%
  x = x + d*.45; %%%%%% additive impulse (long-tailed) noise
subplot(412); plot(d(3001:4000));
down = 1; %%%%%% downsampling
% x = filter(ones(1,down^2)/down^2, 1, x);  % basic lowpass filter
x = x(1:down:end);   % decimation w. or w/o. lp filtering
subplot(413); plot(x(round((3001:4000)/down)));
%  b = ones(1,13); b = b/sum(b); a=1;
%  b = 1; a=[1 -.9];
%  b = 1; a=[1 .9];
%  x = filter(b,a,x);
%  x = medfilt1(x,3);
subplot(414); plot(x(round((3001:4000)/down))); drawnow
%  Fs = 22050; sound(x,Fs);
  Fs = 22050/down; sound(x,Fs);
%  Fs = 22050/down; sound(x*down,Fs);
%  Fs = 22050/down; soundsc(x,Fs);


%%  echo (single)
 x = audioread('dsp01.wav'); Fs = 44100;
% x = audioread('loonytun.wav'); Fs = 11000; 
L = length(x);
% sound(x,Fs);  pause(L/Fs);

path_length = 20; % m
speed = 330; % m/s
delay = path_length / speed * 1000; % ms
dl = round(delay * Fs / 1000); % delay in samples (ms --> samples)
xd = [zeros(dl,1); x(1:end-dl)];
y = .5 * (x + xd);
 sound(y,Fs);
%%% note we could also use y = filter([.5 0 ... 0 .5],1,x);


%%  (partial) left-to-right panning
x = audioread('loonytun.wav'); Fs = 11000;
L = length(x);
xr = x;  xl = x;
mn = 0.3; mx = 1;
xr = xr .* linspace(mn,mx,L)' ./(mn+mx);
xl = xl .* linspace(mx,mn,L)' ./(mn+mx);
sound([xl xr],Fs);


%%  left-vs-right delay  (3.3 ms = 1 m)
x = audioread('loonytun.wav'); % 11000 Hz, 8 bit, 71924 samples
xr = x;  xl = x;
Fs = 11000;
% sound([xl xr],Fs);  pause(L/Fs);

dl = 2 * Fs / 1000; % delay (ms --> samples) % 1 or 2 ms is enough, >10 sounds like reverb
xr = [zeros(dl,1); xr(1:end-dl)];
sound([xl xr],Fs); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   then using Audacity:
% ----- Effect
% echo, reverb
% filter curve eq, graphic eq., bass&treble
% distortion
% noise reduction (requires noise analysis)
%
% ----- Generate
% noise (different freq.content)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
