% E03_4    16mar17   25mar22
% synthesis and "real-time" analysis of sound
% STFT of a chirp

clear all
return

%% synthesis and "real-time" analysis of sound
%%%  NOTE: same-time sound generation and recording may not work on some PCs

Fs = 8000;     % sampling frequency, Hz
dur = .7;       % sound duration, seconds
L = Fs .* dur;  % no. of samples
f_axis = linspace(0,Fs/2,L/2);
rec = audiorecorder(Fs,8,1);
figure(2);

r_all = [];
freq = 2.^[10/12  12/12  8/12  8/12-1  3/12] * 440
for i = 1:5             % 5 short stationary signals
    s = sin(2*pi*freq(i)*dur*[0:L-1]./L);
    s = s.^1; % (odd) power yields harmonic distortion
    sound(s,Fs);
 	recordblocking(rec,dur);
	r = getaudiodata(rec);  % note is a column vector
    r_all = [r_all; r];     % overall non-stationary signal
    R = abs(fft(r));
	subplot(5,1,i); plot(f_axis,R(1:L/2));
    % axis([0 Fs/2 0 40]); 
    drawnow;
	end;
    
R_all = abs(fft(r_all));   % transform of overall non-stationary signal
f_axis_all = linspace(0,Fs/2,5*L/2);
figure(3); plot(f_axis_all,R_all(1:5*L/2));  % sound location in time is unknown
% axis([0 Fs/2 0 40]);

return;

%% same, try shorter duration [without recording]

Fs = 8000;     % sampling frequency, Hz
dur = .8;       % sound duration, seconds
% dur = .2;    % short sounds (.5 and [much] less) produce perceivable 'transition noise' in S_all
L = Fs .* dur;  % no. of samples
f_axis_all = linspace(0,Fs/2,5*L/2);
s_all = [];
freq = 2.^[10/12  12/12  8/12  8/12-1  3/12] * 440  % [Hz]
for i = 1:5             % 5 short stationary signals
    s = sin(2*pi*freq(i)*dur*[0:L-1]./L);
%    s = s + 1*s.^3; % (odd) power yields harmonic distortion
    s_all = [s_all  s];     % overall non-stationary signal
end
sound(s_all,Fs);
S_all = abs(fft(s_all));
figure(1); clf; plot(f_axis_all,S_all(1:5*L/2));  % sound location in time is unknown
figure(2); plot(s_all);


%%  chirp signal: DFT 
clear all;
ome0 = 3e-5*pi;
N = 5000;
x = cos(ome0*[0:N-1].^2.1); % ^2 = linear chirp (try e.g. 1.8, 2.1)
%  W = 200; w = [1./(1+exp(-[-W:W]/100))  ones(1,N-4*W-2)   1./(1+exp(-[W:-1:-W]/100))];
%  x = x .* w;   % windowing the data reduces the Gibbs effect
figure(1); plot(x)
X = abs(fft(x));
figure(2); plot(linspace(0,pi,N/2)/pi, X(1:N/2));   % full DFT

%%        ...and STFT (no overlap, rectangular window, K segments)
K = 5;  % N/K must be an integer; try e.g. 5,25,50
Xstft = [];
figure(1);
for k = 0:K-1,
    xx = x([1:N/K]+k*N/K);  % signal segment
    plot(xx); pause(.9);
    XX = abs(fft(xx))';   % local DFT, (column vector)
    Xstft = [Xstft  XX(1:N/K/2)];  % matrix of DFTs [0...pi]
end
plot(linspace(0,pi,N/K/2)/pi, Xstft);  % per-column plot
figure(2); imshow(flipud(Xstft),[],'InitialMagnification', 200);
      % norm. frequency vs. time  (y-axis ticks are wrong)

%%  STFT with a single Matlab command...
figure(3); spectrogram(x,'yaxis'); % 'yaxis' puts freq. on y axis
figure(4); spectrogram(x,100,'yaxis'); % segment length

return

% help spectrogram;
%    S = spectrogram(X) returns the spectrogram of the signal specified by
%    vector X in the matrix S. By default, X is divided into eight segments
%    with 50% overlap, each segment is windowed with a Hamming window. The
%    number of frequency points used to calculate the discrete Fourier
%    transforms is equal to the maximum of 256 or the next power of two
%    greater than the length of each segment of X.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
 R = audiorecorder;
 record(R);
 pause(2);
 stop(R);
 play(R);
 myR = getaudiodata(R);
 plot(myR);
%%

%%function varargout = GuitarTuner(varargin)
% GUITARTUNER M-file for GuitarTuner.fig
%		
%	Ryan Kinnett, January 2010	

%% AutoTune toy
% https://www.mathworks.com/matlabcentral/fileexchange/26337-autotune-toy

