function spectrogram_FM

	% il segnale oggetto di questo studio ha l'espressione
	%
	% s(t) = cos(2*pi*f0*t + beta * cos(2*pi*fm*t))
	%
	% la banda di questo segnale e' approssimativamente
	% 2*(1+beta)*fm, centrata intorno alla frequenza f0
	%
	ex = input('scegli l''esempio da 0 a 4: ');
	switch ex
		case 0
			fm = 1;
			T_FFT = 1/(16*fm);
		case 1
			fm = 4;
			T_FFT = 1/(4*fm);
		case 2
			fm = 16;
			T_FFT = 1/fm;
		case 3
			fm = 64;
			T_FFT = 4/fm;
		case 4
			fm = 256;
			T_FFT = 16/fm;
	end

	% scegli beta in modo da ottenere la banda richiesta
	B = 2048;
	beta = B/(2*fm)-1;
	
	% durata complessiva del segnale
	T = 2;
	
	% scegli la frequenza di campionamento e la frequenza
	% centrale f0 in modo da soddisfare la condizione di
	% Nyquist e minimizzare le frequenze "inutilizzate"
	fs = 4*B;
	f0 = fs/4;
	fprintf('f0: %g\n',f0);

	% passo di campionamento e numero complessivo di punti
	% della sequenza campionata
	Ts = 1/fs;
	N = T/Ts;
	
	% costruisci la sequenza dei campioni
	t = (0:N-1)*Ts;
	x = cos(2*pi*f0*t + beta * cos(2*pi*fm*t));
	
	% trova N_FFT come potenza di due tale da ottenere la
	% durata richiesta
	N_FFT = 2^ceil(log2(T_FFT/Ts));
	
	% fattore di oversampling per la fft con zero padding
	OSF = 8;
	fprintf('N_FFT: %d\n',N_FFT);

	% vettore delle frequenze discrete
	f = (-N_FFT*OSF/2:N_FFT*OSF/2-1)/(N_FFT*OSF*Ts);
	
	% finestra Gaussiana
	tw = (0:N_FFT-1)-(N_FFT-1)/2;
	gw_truncerr = 1e-2;
	w = exp(log(gw_truncerr)*(tw/tw(end)).^2);
	A_max = sum(abs(w))/2;
	
	% spettrogramma calcolato su N_FFT campioni a passo N_step
	N_step = N_FFT/2;
	NI = ceil((N-N_FFT)/N_step);
	I = zeros(NI,N_FFT*OSF/2);
	
	% offset iniziale
	ofs = 0;
	while ofs+N_FFT <= N
		
		% spettro stimato con la DFT della sequenza finestrata
		a = x(ofs+(1:N_FFT));
		A = fft(a.*w,N_FFT*OSF);

		% passa al modulo e aggiorna l'immagine dello
		% spettrogramma
		tmp = abs(A);
		I = [I(2:end,:);tmp(1:N_FFT*OSF/2)];

		% modulo dello spettro per frequenze positive 
		subplot(2,1,1)
		plot(f,fftshift(tmp));
		set(gca,'XLim',[0 f(end)])
		set(gca,'YLim',[0,A_max])
		ylabel('|A|^2');
		xlabel('f [Hz]');

		% spettrogramma
		subplot(2,1,2)
		imagesc(f(N_FFT*OSF/2+1:end),0:NI-1,I,[0 A_max])
		ylabel('time')
		xlabel('f [Hz]');
		drawnow
		
		% prepara al prossimo loop
		ofs = ofs + N_step;
	end
end
