% E02_1.m             06mar23
%  basic sequence generation and manipulation
clear all;
return

%% "fake-analog" signals vs. discrete sequences
info = audioinfo('dsp01.wav'), pause()
[x_a,Fs]= audioread('dsp01.wav',[1e4+1 4e4]);
x_a = double(x_a) + randn(size(x_a))./100;   % fake-analog signal
figure(1); plot(x_a);
sound(x_a,.6*Fs);
%%
x_a = x_a(10001:11000); % portion of fake-analog signal
figure(2); plot(x_a); 
%%
x_s = x_a(1:10:end); % discrete-time signal (sampling theorem probably violated)
figure(3); stem(x_s); 
figure(4); plot(x_s); pause()
x_qs = round(x_s .* 32); % digital signal (5+1 bit)
figure(5); stem(x_s*32); hold on; stem(x_qs,'r'); hold off
figure(6); stem(x_s.*32 - x_qs);  % quantization error
pause(), close all

%%  Ensemble average
Ns = 100; % no. of samples [sample]
n = 0:Ns-1;
s = n ./ Ns; % signal
 s = n/6 .* exp(-n/15); % another signal
y = 0;
K = 10;
for i = 1:K
    d = randn(size(s))/10;  % noise
    x = s + d;
    y = y + x;
end
y = y/K;

subplot(211);  stem(n,x);
subplot(212);  stem(n,y);

err = y - s;
[std(d).^2,  std(err).^2]


%% same, using matrices
ss = ones(K,1)*s;
dd = randn(size(ss))/10;
xx = ss + dd;
y = mean(xx);

subplot(211);  stem(n,xx(3,:));
subplot(212);  stem(n,y);

%% sample rate alteration (illustrative only)
xa = sin(.07*pi*[0:.01:100]); % fake-analog
xd1 = xa(100:371:end);
xd2 = xa(100:211:end);
subplot(311); plot(xa); axis([0 10000 -1 1])
subplot(312); stem(xd1); axis([0 28 -1 1])
subplot(313); stem(xd2); axis([0 48 -1 1])
% we could generate xd2 starting from xd1 and viceversa
% without actually resampling xa 
% (in practical cases, different results would be obtained)

%%  downsampling
Ns = 100; % no. of samples [sample]
n = 0:Ns-1;
ome0 = 0.1*pi;
x = sin(ome0*n);
for k = 1:15,
y = x(1:k:end);
Y = abs(fft(y));
subplot(211);stem(y);

subplot(212); stem(Y);   % note aliasing for 'large' k
pause();
end

%%  upsampling (Matlab trick)
xx = [3 6 3 2 7];
xu = [xx; zeros(size(xx)); zeros(size(xx))]
xu = xu(:); xu = xu'
subplot(211); stem(xx);
subplot(212); stem(xu);


%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BASIC SEQUENCES
%%  complex exponential
Ns = 100; % no. of samples [sample]
n = 0:Ns-1;
ome0 = .25*pi;      % ang.freq.  [rad/sample]  [0, pi)     
sigma = -0.1;         % amplitude variation speed  (-S, S)      [+-.1]

c = sigma + j*(ome0);
x = exp(c*n);
% figure(1);
subplot(211);  stem(n,real(x));
xlabel('n'); ylabel('Amplitude'); title('Real');
subplot(212);  stem(n,imag(x));
xlabel('n'); ylabel('Amplitude'); title('Imag');

%%  real sinusoid: true period rarely is 2pi/omega
Ns = 111; % try Ns = integer no. of periods or not
n = 0:Ns-1;
ome0 = .1*pi;  % period = 2pi/ome0 if 2pi/ome0 integer
               % otherwise period = r 2pi/ome0
phi = pi/4; % phase [rad] [-pi pi]
x = sin([ome0*n+phi]);
subplot(211); stem(n,x);
subplot(212); stem(abs(fft(x)));

% note same holds for ome0 = .1*pi+(2*k*pi);
% digital frequency is periodical

%%  real sinusoid: digital frequency is symmetrical wrt pi
Ns = 1000;
n = 0:1:Ns-1;
for ome0 = [.0:.11:10]*pi; % ome0 monot. increases
    x = sin(ome0*n);
    subplot(211); plot(x(1:50));
    subplot(212); plot(abs(fft(x))); axis([0 Ns 0 600]);
    Fs = 4000; sound(x(1:200),Fs); 
    pause(.2)
end
hold off;

