rayleigh_fad.c

ofdm的完整系统模型,包含信道参数,多径模型,doppler频移等都可以自由修改!是您做仿真的有力帮助.c语言运行速度快!

#include "cgs.h"

long	cgs_noiseseed;

void initial_flat(fad_params)
FadParams *fad_params;
{
	double	fr;
	double	*cp, sum, half = 0.5 * (fad_params->P_filter_len-1);
	int	i;
	
	if (strcmp(fad_params->P_interpolate, "yes") == 0 && fad_params->P_int_factor > 1)
		fr = PI2 * fad_params->P_fd * fad_params->P_int_factor / fad_params->P_s_freq;
	else
		fr = PI2 * fad_params->P_fd / fad_params->P_s_freq;
		
	switch(fad_params->P_window)
	{
		case  BLACKMAN: 
			fad_params->S_scoef = DvGenerate(OvBLACKMAN, fad_params->P_filter_len, NULL);
			break;
		case BARTLETT : 
			fad_params->S_scoef = DvGenerate(OvBARTLETT, fad_params->P_filter_len, NULL);
			break;
		case HAMMING  : 
			fad_params->S_scoef = DvGenerate(OvHAMMING, fad_params->P_filter_len, NULL);
			break;
		case HANNING : 
			fad_params->S_scoef = DvGenerate(OvHANNING, fad_params->P_filter_len, NULL);
			break;
		case RECTANGLAR : 
			fad_params->S_scoef = OvAlloc(fad_params->P_filter_len, NULL);
			Dvfill(1.0, fad_params->S_scoef);
			break;
	}
	
	//Generate the impulse response
	cp = (double *) OvGetStart(fad_params->S_scoef);
	#ifdef U_SPECTRUM
	{
		for (i = 0; i < fad_params->P_filter_len; i++) {
			sum = fabs(fr * (i - half));
			if (sum == 0.0) {
				*cp++ *= J25(1.0e-12) * 1000.0;
			}
			else {
				*cp++ *= J25(sum) / pow(sum, 0.25);
			}
		}		
	}
	#else
	{
		double magi = sqrt(2.0* fad_params->P_fd);
		for (i = 0; i < fad_params->P_filter_len; i++) {
			sum = fabs(fr * (i - half));
			if (sum == 0.0) {
				*cp++ *= magi;
			}
			else {
				*cp++ *= magi * sin(sum) / sum ;
			}
		}
	}
	#endif
	
	
	if(strcmp(fad_params->P_norm, "yes") == 0) {
   		Dvdot(fad_params->S_scoef, fad_params->S_scoef, &sum);
		Dvscale(sqrt(0.5 * fad_params->P_norm_power/sum), fad_params->S_scoef, fad_params->S_scoef);
	}
	else{
		Dvcsum(fad_params->S_scoef, &sum);
		Dvscale(sqrt(0.5) * fad_params->P_gain/sum, fad_params->S_scoef, fad_params->S_scoef);
	}


	fad_params->S_stater = OvAlloc(fad_params->P_filter_len, NULL);
	fad_params->S_statei = OvAlloc(fad_params->P_filter_len, NULL);
	
	if(strcmp(fad_params->P_interpolate, "yes") == 0 && fad_params->P_int_factor > 1) {
		fad_params->S_intrp = OvAlloc(fad_params->P_nipts * fad_params->P_int_factor, NULL);
		filtGenWSincLpf(fad_params->S_intrp, PI / fad_params->P_int_factor, (double) fad_params->P_int_factor);

		OvSetStepSize(fad_params->S_intrp, fad_params->P_int_factor);
		OvSetLength(fad_params->S_intrp, fad_params->P_nipts);
		fad_params->S_istater = OvAlloc(fad_params->P_nipts, NULL);
		fad_params->S_istatei = OvAlloc(fad_params->P_nipts, NULL);
		fad_params->S_cnt = 0;
	}
	
/*
 *  Fill up the state vectors
 */
	{
		double	x;
		double	*sr = (double *) OvGetStart(fad_params->S_stater);
		double	*si = (double *) OvGetStart(fad_params->S_statei);
		int	i, n = OvGetLength(fad_params->S_stater);

		for (i = 0; i < n; i++) {
			*sr++ = noiseWGN(&cgs_noiseseed);
			*si++ = noiseWGN(&cgs_noiseseed);
		}
		
		if (strcmp(fad_params->P_interpolate, "yes") == 0 && fad_params->P_int_factor > 1) {
			sr = ((double *) OvGetEnd(fad_params->S_istater)) - 1;
			si = ((double *) OvGetEnd(fad_params->S_istatei)) - 1;
			for (i = 0; i < fad_params->P_nipts; i++) {
				x = noiseWGN(&cgs_noiseseed);
				*sr-- = Dvfir(x, fad_params->S_stater, fad_params->S_scoef);
				x = noiseWGN(&cgs_noiseseed);
				*si-- = Dvfir(x, fad_params->S_statei, fad_params->S_scoef);
			}
		}
	}

}


Complex rayleigh_flat(cpx_in, fad_params)
Complex cpx_in;
FadParams *fad_params;
{

	Complex	w, fad_out;
	double	v;
	
	if(strcmp(fad_params->P_interpolate, "yes") == 0 && fad_params->P_int_factor > 1) {
		if (fad_params->S_cnt-- < 2) {
			fad_params->S_cnt = fad_params->P_int_factor;
			v = noiseWGN(&cgs_noiseseed);
			w.real = Dvfir(v, fad_params->S_stater, fad_params->S_scoef);
			v = noiseWGN(&cgs_noiseseed);
			w.imag = Dvfir(v, fad_params->S_statei, fad_params->S_scoef);
			OvSetVirtStart(fad_params->S_intrp, OvGetStart(fad_params->S_intrp));
			w.real = Dvfir2(w.real, fad_params->S_istater, fad_params->S_intrp);
			w.imag = Dvfir2(w.imag, fad_params->S_istatei, fad_params->S_intrp);
		}
		else {
			OvSetVirtStart(fad_params->S_intrp, OvGetVirtStart(fad_params->S_intrp) + sizeof(double));
			Dvdot(fad_params->S_istater, fad_params->S_intrp, &w.real);
			Dvdot(fad_params->S_istatei, fad_params->S_intrp, &w.imag);
		}             
	}
	else {
		v = noiseWGN(&cgs_noiseseed);
		w.real = Dvfir(v, fad_params->S_stater, fad_params->S_scoef);
		v = noiseWGN(&cgs_noiseseed);
		w.imag = Dvfir(v, fad_params->S_statei, fad_params->S_scoef);	
	}
		
		
	fad_out.real = w.real * cpx_in.real - w.imag * cpx_in.imag;
	fad_out.imag = w.real * cpx_in.imag + w.imag * cpx_in.real;
	
	return fad_out;

}

void initial_sel(ray_params)
RayParams *ray_params;
{
	int i, nv_delay[12];
	double f_doppler, absv_pow[12], total_pow;
	long interpolater_len;
	
	/* initial noise seed for the same channel of each loop */
	cgs_noiseseed = 1550;

	total_pow = 0;
	
	// delay and power of each path
	#ifdef CHNANNEL_COST207
	{
		double factor_a;
		factor_a = 3 * log(10.0) / MAX_DELAY;
		for (i = 0; i < NUM_MULTIPATH; i++) {
			nv_delay[i] = (int) floor(i * MAX_DELAY * CH_SAMPLE_RATE / (NUM_MULTIPATH-1) + 0.5);
			//absv_pow[i] = factor_a * exp(-factor_a*nv_delay[i]/CH_SAMPLE_RATE) / (1-exp(-factor_a*MAX_DELAY));
			absv_pow[i] = exp(-factor_a*nv_delay[i]/CH_SAMPLE_RATE);	// power relative to first path
			total_pow = total_pow + absv_pow[i];
		}
		for (i = NUM_MULTIPATH; i < 12; i++) {
			nv_delay[i] = 0;
			absv_pow[i] = 0;
		}
	}
	#else
	{
		for (i = 0; i < NUM_MULTIPATH; i++) {
			nv_delay[i] = (int) floor(vec_delay_t[i] * CH_SAMPLE_RATE + 0.5);
			absv_pow[i] = pow(10.0, vec_power_db[i]/10.0);
			total_pow = total_pow + absv_pow[i];
		}
		for (i = NUM_MULTIPATH; i < 12; i++) {
			nv_delay[i] = 0;
			absv_pow[i] = 0;
		}
	}
	#endif
	
	// fading channel normalized gain
	ray_params->Ch_Norm_Factor = 1.0 / sqrt(total_pow);
	
	// doppler frequency
	#ifdef FD_CACULATE
		f_doppler = (CH_CARRIER_FREQ * CH_MOBILE_SPEED) / (3.6 * 3.0e8);
	#else
		f_doppler = CH_DOPPLER_FREQ;
	#endif
	
	// interpolater length
	#ifdef FD_NORMALIZE
		interpolater_len = (long) floor(FD_NORMALIZE*CH_SAMPLE_RATE/f_doppler + 0.5);
	#else
		interpolater_len = CH_INTERPOLATE_LEN;
	#endif		
	
	// initial each filter's paramter
	for (i = 0; i < NUM_MULTIPATH; i++) {
		ray_params->fad_params[i].delay_n = nv_delay[i];
		ray_params->fad_params[i].P_s_freq = CH_SAMPLE_RATE;
		ray_params->fad_params[i].P_fd = f_doppler;
		ray_params->fad_params[i].P_int_factor = interpolater_len;
		ray_params->fad_params[i].P_filter_len = CH_FILTER_LEN;
		ray_params->fad_params[i].P_nipts = CH_INTERPOLATE_NUM;
		ray_params->fad_params[i].P_gain = sqrt(absv_pow[i]);
		ray_params->fad_params[i].P_norm_power = absv_pow[i];
		ray_params->fad_params[i].P_window = CH_WINDOW;
		strcpy(ray_params->fad_params[i].P_norm, CH_IF_NORMALIZE);
		strcpy(ray_params->fad_params[i].P_interpolate, CH_IF_INTERPOLATE);
	}

	// initial each filter's state 
	for (i = 0; i < NUM_MULTIPATH; i++) {
		initial_flat(&(ray_params->fad_params[i]));
	}
	
	// initial circular buffer
	ray_params->Data_Buff = OvAlloc((nv_delay[NUM_MULTIPATH-1] + 1)*2, NULL);
	
	// logfile
	#ifdef LOG_SCOEF_FILE
	{
		FILE *fp_scoef;
		//char *log_scoef = LOG_SCOEF_FILE;
		double *coef_p = (double *) OvGetStart(ray_params->fad_params[0].S_scoef);
		
		fp_scoef = fopen(LOG_SCOEF_FILE, "w");
		for (i = 0; i < OvGetLength(ray_params->fad_params[0].S_scoef); i++) {
			fprintf(fp_scoef, "%f\n", *coef_p);
			coef_p++;
		}
		fclose(fp_scoef);
	}
	#endif
	
}


Complex rayleigh_sel(cpx_in, ray_params)
Complex		cpx_in;
RayParams	*ray_params;
{
	int i;
	Complex sum, filt_out;
	
	register Complex *sp = (Complex *) ray_params->Data_Buff->curr_loc;
	register Complex *se = (Complex *) ray_params->Data_Buff->virtend;
	register Complex *cur_idx = (Complex *) ray_params->Data_Buff->curr_loc;
	
 	// update the circular buffer pointer.
	if (sp-- <= (Complex *) ray_params->Data_Buff->virtstart) {
		sp = se - 1;
	}
	*sp = cpx_in;
	ray_params->Data_Buff->curr_loc = (char *) sp;
	
	// pass to channel
	sum.real = 0;
	sum.imag = 0;
	
	for (i = 0; i < NUM_MULTIPATH; i++) {
		cur_idx = sp + ray_params->fad_params[i].delay_n;
		if (cur_idx >= se) cur_idx = ((Complex *) ray_params->Data_Buff->virtstart) + (cur_idx - se);
		filt_out = rayleigh_flat(*cur_idx, &(ray_params->fad_params[i]));
		sum = c_add(sum, filt_out);
	}
	sum.real = sum.real * ray_params->Ch_Norm_Factor;
	sum.imag = sum.imag * ray_params->Ch_Norm_Factor;

	return sum;
}


/* free the memory of the rayleigh channel */
void free_rayleigh_sel(ray_params)
RayParams	*ray_params;
{
	int i;

	/* free the multipath buffer */
	if (NULL != ray_params->Data_Buff->start) {
		free(ray_params->Data_Buff->start);
		ray_params->Data_Buff->start = NULL;
		ray_params->Data_Buff->virtstart = NULL;
	}

	/* free the fadding channel buffer */
	for(i = 0; i < NUM_MULTIPATH; i++) {
		
		if (NULL != ray_params->fad_params[i].S_scoef->start) {
			free(ray_params->fad_params[i].S_scoef->start);
			ray_params->fad_params[i].S_scoef->start = NULL;
			ray_params->fad_params[i].S_scoef->virtstart = NULL;
		}		
		if (NULL != ray_params->fad_params[i].S_stater->start) {
			free(ray_params->fad_params[i].S_stater->start);
			ray_params->fad_params[i].S_stater->start = NULL;
			ray_params->fad_params[i].S_stater->virtstart = NULL;
		}
		if (NULL != ray_params->fad_params[i].S_statei->start) {
			free(ray_params->fad_params[i].S_statei->start);
			ray_params->fad_params[i].S_statei->start = NULL;
			ray_params->fad_params[i].S_statei->virtstart = NULL;
		}
		if (NULL != ray_params->fad_params[i].S_istater->start) {
			free(ray_params->fad_params[i].S_istater->start);
			ray_params->fad_params[i].S_istater->start = NULL;
			ray_params->fad_params[i].S_istater->virtstart = NULL;
		}
		if (NULL != ray_params->fad_params[i].S_istatei->start) {
			free(ray_params->fad_params[i].S_istatei->start);
			ray_params->fad_params[i].S_istatei->start = NULL;
			ray_params->fad_params[i].S_istatei->virtstart = NULL;
		}
		if (NULL != ray_params->fad_params[i].S_intrp->start) {
			free(ray_params->fad_params[i].S_intrp->start);
			ray_params->fad_params[i].S_intrp->start = NULL;
			ray_params->fad_params[i].S_intrp->virtstart = NULL;
		}
	}

}

/* initial AWGN channel */
void initial_gauss(gauss_params, chn_delay, signal_pow, chn_snr)
GauParams *gauss_params;
double chn_delay;
double signal_pow;
double chn_snr;
{
	
	gauss_params->Ch_Delay = (int) floor(chn_delay * CH_SAMPLE_RATE + 0.5);
	gauss_params->Norm_Var = sqrt(0.5 * signal_pow / pow(10.0, chn_snr/10.0));
	gauss_params->Delay_Buff = OvAlloc(2*(gauss_params->Ch_Delay+1), NULL);
	Dvfill(0.0, gauss_params->Delay_Buff);
}

/* AWGN Channel */
Complex awgn_channel(cpx_in, gauss_params)
Complex	cpx_in;
GauParams *gauss_params;
{
	double noise;
	Complex cpx_out, cpx_tmp;

	register Complex *sp = (Complex *) gauss_params->Delay_Buff->curr_loc;
	register Complex *se = (Complex *) gauss_params->Delay_Buff->virtend;
	register Complex *cur_idx = (Complex *) gauss_params->Delay_Buff->curr_loc;

 	// update the circular buffer pointer.
	if (sp-- <= (Complex *) gauss_params->Delay_Buff->virtstart) {
		sp = se - 1;
	}
	*sp = cpx_in;
	gauss_params->Delay_Buff->curr_loc = (char *) sp;

	/* read data from delay_buffer */
	cur_idx = sp + gauss_params->Ch_Delay;
	if (cur_idx >= se) cur_idx = ((Complex *) gauss_params->Delay_Buff->virtstart) + (cur_idx - se);
	cpx_tmp = *cur_idx;
	
	noise = gauss_params->Norm_Var * noiseWGN(&cgs_noiseseed);
	cpx_out.real = cpx_tmp.real + noise;
	noise = gauss_params->Norm_Var * noiseWGN(&cgs_noiseseed);
	cpx_out.imag = cpx_tmp.imag + noise;
	
	return cpx_out;
}

/* free the memory of the AWGN channel */
void free_gauss(gauss_params)
GauParams *gauss_params;
{
	if (NULL != gauss_params->Delay_Buff->start) {
		free(gauss_params->Delay_Buff->start);
		gauss_params->Delay_Buff->start = NULL;
		gauss_params->Delay_Buff->virtstart = NULL;
	}
}


/* initial carrier offset channel */
void initial_carrier_offset(freq_params, freq_shift, phase_ini)
FreqParams *freq_params;
double freq_shift;
double phase_ini;
{
	freq_params->Phase_Reg = phase_ini;
	freq_params->Base_Phase = PI2*freq_shift/SUBCARRIERS;			// relative to subcarrier
//	freq_params->Base_Phase = PI2*freq_shift/CH_SAMPLE_RATE;
}


/* Frequency Offset */
Complex carrier_offset(cpx_in, freq_params)
Complex	cpx_in;
FreqParams *freq_params;
{
	Complex cpx_out;
	Complex phase_rotate;

	phase_rotate.real = cos(freq_params->Phase_Reg);
	phase_rotate.imag = sin(freq_params->Phase_Reg);

	cpx_out = c_mul(cpx_in, phase_rotate);

	freq_params->Phase_Reg = freq_params->Phase_Reg + freq_params->Base_Phase;
	if (fabs(freq_params->Phase_Reg) >= PI2) {
		freq_params->Phase_Reg = (freq_params->Phase_Reg > 0)? freq_params->Phase_Reg-PI2 : freq_params->Phase_Reg+PI2;
	}
	
	return cpx_out;
}