clear all;
k				= 1.381e-23;	% Joule/K 	- Boltzmann's Constant
T0				= 290;			% K			- temperature
% Encoding Windows-1251

% 1.    Функциональная схема
% 1010          -   Схема защиты входных  каскадов
% 1011          -   Преселектор 
% 1012          -   Управляемый МШУ (вкл/выкл.)
% 1013          -   Фильтр основной селекции 
% 1014          -   Управляемый усилитель петли АРУ
% 1015          -   Квадратурный смеситель
% 1021, 1022    -   Фильтры низкой частоты.
% 1031, 1032    -   Низкочастотные усилители
% 1041, 1042    -   Фильтры низкой частоты АЦП
% 1051, 1052    -   АЦП

% 2.    Системная диаграмма радиотракта для анализа мощности шумов и чувствительности

		Rrf_inp_ohm			= 50;		% Ohm	- 	входной импеданс.

% 2011  -   Преселектор 
        BWrf_ekv_prf_hz   	= 20.0e6;   % Hz	-	Эквивалентная полоса пропускания преселектора 
        Lrf_max_prf_db      =  1.0;     % dB 	- 	Максимальный уровень затухания преселектора в рабочей полосе РПУ (insertion loss)
        
% 2012  -   МШУ 
% (для примера взяты параметры Agilent MGA-71543 с учетом реализационных потерь)
        Grf_lna_db     	    = 16.0;     % dB 	- 	Коэффициент усиления МШУ 
        NFrf_lna_db         = 1.0;      % dB		- 	Коэффициент шума МШУ 
        
% 2013  -   Фильтр основной селекции 
        BWrf_ekv_fms_hz   	=  6.0e6;   % Hz	-	Полоса пропускания фильтра основной селекции 
        Lrf_max_fms_db      =  4.0;    % dB 	- 	Максимальный уровень затухания фильтра основной селекции в полосе пропускания (insertion loss)
        
% 2014  -   Усилитель высокой частоты 
% (для примера взяты параметры Agilent MGA-71543 с учетом реализационных потерь)
        Grf_amp_db      	=  16.0;    % dB 	- 	Коэффициент усиления УВЧ 
        NFrf_amp_db     	=  1.0;     % dB    - 	Коэффициент шума УВЧ
		
% 2015  -   Квадратурный смеситель
% (для примера взяты параметры Analog Devices ADL5387)
        Grf_mix_db      	=   4.5;    	% dB 	- 	Коэффициент усиления квадратурного смесителя 
        NFrf_mix_db     	=  15.0;    	% dB    - 	Коэффициент шума квадратурного смесителя
		IP1dBrf_mix_dbw		=  13.0 - 30.0; % dBW	- 	Input P1dB (IP1dB)

% 3.    Системная диаграмма цифрового приемника для анализа чувствительности и максимально допустимого уровня блокирующего сигнала
% 3011	-	Радиотракт
% Коэффициент шума радиотракта (db)
        NFrf_sys_db = pow2db( ( db2pow(  Lrf_max_prf_db )       ) + ...                   
                              ( db2pow( NFrf_lna_db     ) + 1 ) / ( db2pow( -Lrf_max_prf_db )) + ...
                              ( db2pow(  Lrf_max_fms_db ) + 1 ) / ( db2pow( -Lrf_max_prf_db ) * db2pow(  Grf_lna_db )) + ...
                              ( db2pow( NFrf_amp_db     ) + 1 ) / ( db2pow( -Lrf_max_prf_db ) * db2pow(  Grf_lna_db ) * db2pow( -Lrf_max_fms_db )) + ...
							  ( db2pow( NFrf_mix_db     ) + 1 ) / ( db2pow( -Lrf_max_prf_db ) * db2pow(  Grf_lna_db ) * db2pow( -Lrf_max_fms_db ) * db2pow(  Grf_amp_db )) ...
                            );

% Коэффициент усиления радиотракта (dB)
		Grf_sys_db	= ( Grf_lna_db + Grf_amp_db + Grf_mix_db ) - ( Lrf_max_prf_db + Lrf_max_fms_db );
		
% Полоса сигнала на выходе радиотракта (Hz)
		BWrf_sys_hz = BWrf_ekv_fms_hz;

% 3012	-	Baseband LPF
		Lbb_lpf_db		= 9;		% dB			-	Максимальный уровень затухания Baseband LPF в полосе пропускания (insertion loss)
		
% 3013	-	Baseband УНЧ (LTC6400-14)
		Gbb_opa_db		= 0;		% dB			-	Коэффициент усиления УНЧ
		NFbb_opa_db     = 0;      % dB			-	Коэффициент шума УНЧ	

% Эквивалентный коэффициент усиления УНЧ, состоящего из фильтра и ОУ (dB)
		Gbb_lfa_db		= Gbb_opa_db - Lbb_lpf_db;
		
% Коэффициент шума системы от входа до выхода УНЧ (dB) - практически не поменялся. 
		NFbb_sys_db     = pow2db( ( db2pow( NFrf_sys_db         )     ) + ...                   
                                  ( db2pow( NFbb_opa_db         ) + 1 ) / ( db2pow( Grf_sys_db )) ...
					    		);
							
% Коэффициент усиления системы от входа до выхода УНЧ (dB)
		Gbb_sys_db	    = Grf_sys_db + Gbb_lfa_db; 

% Полоса сигнала  в основной полосе (Hz) - не поменялась
		BWbb_sys_hz		= BWrf_sys_hz;
		
% Мощность шума на выходе УНЧ в основной полосе (dBW)		
		PNbb_out_dbw	= pow2db( k * T0 * BWrf_sys_hz ) + NFbb_sys_db + Gbb_sys_db;	

% 3014	-	АЦП 
% (для примера взяты параметры 1 канала Linear Technology LTC2271)
		FSadc_hz		= 20.0e6;	% Hz	-	Sampling rate		
		SNadc_fs_db		= 84;		% dB	-	SNR 
		NBadc_fs_bits	= 16;		% bits 	-	Full scale bits
		Vadc_fs_v		= 2;		% V		-	Full scale voltage
		Radc_inp_ohm	= 1000;		% Ohm	-	Input ADC resistance	

% Мощность синусоидального сигнала при размахе на полную шкалу для пары АЦП (dBW)			
		PFSadc_inp_dbw	= pow2db( 2.0 * (( Vadc_fs_v * 0.5 * sqrt( 0.5 )) ^ 2 ) / Radc_inp_ohm );  
		
% Мощность собственного шума пары АЦП (dBW)			
		PNadc_snr_dbw	= PFSadc_inp_dbw - SNadc_fs_db;
		
% Мощность шума квантования пары АЦП (dBW)
		PNadc_qan_dbw	= PFSadc_inp_dbw - ( NBadc_fs_bits * mag2db( 2 ) + mag2db( sqrt( 6 ) / 2 )); % adding correction factor for sinusoidal signal			
					
% 3015	-	Эквивалентная схема обработки сигнала с ЧМ
        SNfm_min_db		= 12.0;		% dB  		- минимальное отношение сигнал-шум для приема голосовых передач в режиме ЧМ
		BWfm_max_hz		= 25.0e3;	% Hz 		- максимальная ширина полосы канала, предназначенного для передачи речи		
		BWfm_min_hz		= 6.25e3;	% Hz 		- минимальная   ширина полосы канала, предназначенного для передачи речи

% Чувствительность радиотракта (dBW) при идеальных АЦП и цифровой обработке канала передачи (минимальная мощность в полосе канала передачи)	
		Pfm_min_dbw		= pow2db( k * T0 * BWfm_max_hz ) + NFbb_sys_db + SNfm_min_db;		
		
% Мощность сигнала на выходе унч при синусоидальном входном сигнале, по уровню равным уровню чувствительности (dBW)
		Pfm_min_bb_sys_dbw 	= Pfm_min_dbw + Gbb_sys_db;		

% Мощность сигнала на входе смесителя. Такой сигнал после прохождения всех последующих каскадов 
% на выход АЦП даст размах полной шкалы (dBW)
		PFSmix_inp_dbw = PFSadc_inp_dbw - Gbb_lfa_db - Grf_mix_db;		

% Запас линейного участка передаточной характеристики смесителя (по входу)
		deltaPmix_inp_lin = IP1dBrf_mix_dbw - PFSmix_inp_dbw; 
		

% Разрешение БПФ равно минимальной полосе канала 
Nfft = 2 * (( FSadc_hz / 2 ) / BWfm_min_hz ); 
Nsmp = Nfft;

tmp_fft_buf = zeros( 1, Nfft );
tmp_acc_buf = zeros( 1, Nfft );
tmp_smp_buf = zeros( 1, Nsmp );

max_acc = 30;

for acc = 1:max_acc
% Принимаемый сигнал 1 
% на входе приемного тракта - сигнал на уровне чувствительности приемного тракта
% на входе АЦП - сигнал прошедший все каскады преобразования и усиления

  PS1 = db2pow( Pfm_min_bb_sys_dbw );
  WS1 = 25.0e3;
  FS1 =  1.0e6;
  
  Fstart = FS1;
  Fstop  = Fstart + WS1 - BWfm_min_hz;
  Pstep  = PS1 / ( WS1 / BWfm_min_hz ); 
  PS1_smp_buf = zeros( 1, Nsmp );
  
  for f = Fstart:BWfm_min_hz:Fstop
      phi_acc = 2.0 * pi * rand( 1 ); % random phase
      phi_stp =       pi * ( f / ( FSadc_hz / 2 ));
      
      for k = 1:Nsmp
          PS1_smp_buf( k ) = PS1_smp_buf( k ) + sqrt( Pstep ) * exp( j * phi_acc );
          phi_acc = phi_acc + phi_stp;
          
          if( phi_acc > (  2.0 * pi ))
              phi_acc = phi_acc - 2.0 * pi;
          else
              if( phi_acc < ( -2.0 * pi ))
                  phi_acc = phi_acc + 2.0 * pi;
              end
          end
      end
  end

% Принимаемый сигнал 2 
% на входе приемного тракта - сигнал по уровню близок блокировке приемного тракта
% на входе АЦП - сигнал прошедший все каскады преобразования и усиления

  PS2 = db2pow( PFSadc_inp_dbw - 1.0 ); % -1 dB back off 
  WS2 = 25.0e3;
  FS2 = -2.0e6;
  
  Fstart = FS2;
  Fstop  = Fstart + WS2 - BWfm_min_hz;
  Pstep  = PS2; 
  PS2_smp_buf = zeros( 1, Nsmp );
  
  for f = Fstart:BWfm_min_hz:Fstop
      phi_acc = 2.0 * pi * rand( 1 ); % random phase
      phi_stp =       pi * ( f / ( FSadc_hz / 2 ));
      
      for k = 1:Nsmp
          PS2_smp_buf( k ) = PS2_smp_buf( k ) + sqrt( Pstep ) * exp( j * phi_acc );
          phi_acc = phi_acc + phi_stp;
          
          if( phi_acc > (  2.0 * pi ))
              phi_acc = phi_acc - 2.0 * pi;
          else
              if( phi_acc < ( -2.0 * pi ))
                  phi_acc = phi_acc + 2.0 * pi;
              end
          end
      end
  end

% Шум на выходе унч основной полосы
  PN1 = db2pow( PNbb_out_dbw );
  WN1 = BWbb_sys_hz;
  
  Pfull_bw = PN1 * ( FSadc_hz / WN1 );
  
  PN1_smp_buf = sqrt( 0.5 * Pfull_bw ) * complex( randn( 1, Nsmp ), randn( 1, Nsmp ));
  tmp_fft_buf = fftshift( fft( PN1_smp_buf ));
  tmp_msk_buf = zeros( 1, Nfft );
  tmp_msk_buf((( Nfft / 2 ) - (( WN1 / FSadc_hz ) * ( Nfft / 2 )) + 1 ) : (( Nfft / 2 ) + (( WN1 / FSadc_hz ) * ( Nfft / 2 )))) = ... 
              ones( 1, (( WN1 / FSadc_hz ) * Nfft ));
  
  tmp_fft_buf = tmp_fft_buf .* tmp_msk_buf;
  PN1_smp_buf = ifft( fftshift( tmp_fft_buf ));

% Шум аналоговой части АЦП
  PN2 = db2pow( PNadc_snr_dbw ) - db2pow( PNadc_qan_dbw );
  %PN2 = db2pow( PNadc_snr_dbw );
  
  Pfull_bw = PN2;
  
  PN2_smp_buf = sqrt( 0.5 * Pfull_bw ) * complex( randn( 1, Nsmp ), randn( 1, Nsmp ));

% Ограничитель и квантователь АЦП
  QAN_smp_buf = PS1_smp_buf + PS2_smp_buf + PN1_smp_buf + PN2_smp_buf;
  
  QAN_delta = Vadc_fs_v / ( 2 ^ NBadc_fs_bits );
  
  QAN_smp_buf = round( QAN_smp_buf ./ QAN_delta ) .* QAN_delta;
  
  QAN_smp_buf_re = real( QAN_smp_buf );
  QAN_smp_buf_re( find( QAN_smp_buf_re > (  Vadc_fs_v / 2.0 ))) =  Vadc_fs_v / 2.0;
  QAN_smp_buf_re( find( QAN_smp_buf_re < ( -Vadc_fs_v / 2.0 ))) = -Vadc_fs_v / 2.0;
  
  QAN_smp_buf_im = imag( QAN_smp_buf );
  QAN_smp_buf_im( find( QAN_smp_buf_im > (  Vadc_fs_v / 2.0 ))) =  Vadc_fs_v / 2.0;
  QAN_smp_buf_im( find( QAN_smp_buf_im < ( -Vadc_fs_v / 2.0 ))) = -Vadc_fs_v / 2.0;
  
  QAN_smp_buf = complex( QAN_smp_buf_re, QAN_smp_buf_im );

% Спектр после ADC

  tmp_smp_buf = QAN_smp_buf;
  %tmp_smp_buf = PS1_smp_buf + PS2_smp_buf + PN1_smp_buf + PN2_smp_buf;
  
  tmp_fft_buf = fft( tmp_smp_buf ) / Nfft;
  tmp_acc_buf = tmp_acc_buf + ( tmp_fft_buf .* conj( tmp_fft_buf ));
end  

tmp_acc_buf = tmp_acc_buf ./ max_acc;

f = linspace(( -FSadc_hz / 2 ) + BWfm_min_hz, FSadc_hz / 2, Nfft );

plot( f, pow2db( fftshift( tmp_acc_buf ))); 
xlim( [( -FSadc_hz / 2 ), ( FSadc_hz / 2 )] ); 
ylim( [-150.0, -20.0] ); 
title('Power Spectrum')
xlabel('Frequency (Hz)')
ylabel('P(f) dBW')
drawnow;