DSPRelated.com
Code
   Code Snippets    Submitted by Amir Shafiq

Passive direction finding

● December 15, 2010 ● Coded in Matlab 
clc
clear all
close all

fs=20e7;
T=1/fs;

t=0:T:0.0000002;
f=10e6;
scanagle_deg=60;
step_deg=6;

theta=-deg2rad(scanagle_deg):deg2rad(step_deg):deg2rad(scanagle_deg);

for i=1:length(theta);

ant1=2*sin(2*pi*f*t);
ant2=2*sin(2*pi*f*t+theta(i));

[my,mx]=max(ant1);

sum_ant(i)=ant1(mx)+ant2(mx);
diff_ant(i)=ant1(mx)-ant2(mx);
ratio_ant(i)=diff_ant(i)/sum_ant(i);

% if diff_ant(i)==0
%     diff_ant(i)=0.1;
% end
% ratio1_ant(i)=sum_ant(i)/diff_ant(i);

end

% subplot(311)
% plot(t,ant1,t,ant2)
subplot(211)
plot(rad2deg(theta),sum_ant,rad2deg(theta),diff_ant)
subplot(212)
plot(rad2deg(theta),ratio_ant)

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

sim_ant1=2*sin(2*pi*f*t);
sim_ant2=2*sin(2*pi*f*t+deg2rad(30));

[smy,smx]=max(sim_ant1);

%%%%also take for same sample value of data

sim_sum_ant=sim_ant1(mx)+sim_ant2(mx);
sim_diff_ant=sim_ant1(mx)-sim_ant2(mx);
sim_ratio_ant=sim_diff_ant/sim_sum_ant

passive Sum difference method for finding bearing

● December 15, 2010 ● Coded in Matlab 
clc
clear all
close all

fs=20e7;
T=1/fs;

t=0:T:0.0000002;
f=10e6;
scanagle_deg=60; % Range covered by seeker antenna
step_deg=6;

% for ploting original arrays of sum, diff and ratio
theta=-deg2rad(scanagle_deg):deg2rad(step_deg):deg2rad(scanagle_deg);

for i=1:length(theta);

ant1=1*sin(2*pi*f*t);               %Antenna 1
ant2=1*sin(2*pi*f*t+theta(i));      %Antenna 2

[my,mx]=max(ant1);

sum_ant(i)=ant1(mx)+ant2(mx);               %Sum of antennas
diff_ant(i)=ant1(mx)-ant2(mx);              %diff of antennas
ratio_ant(i)=diff_ant(i)/sum_ant(i);        %ratio

end

% subplot(311)
% plot(t,ant1,t,ant2)
% subplot(211)
% plot(rad2deg(theta),sum_ant,rad2deg(theta),diff_ant)
% subplot(212)
% plot(rad2deg(theta),ratio_ant)

% 
% figure
% subplot(211)
% plot(rad2deg(theta),sumn_ant,rad2deg(theta),diffn_ant)
% subplot(212)
% plot(rad2deg(theta),ration_ant)

%%%%%%%%%%%%%%%%Recieved signal at antenna for DOA%%%%%%%%%%%%%%%%%%%%%%%%%%

A_ant1=2;                   % Amplitude of antenna 1
A_ant2=3;                   % Amplitude of antenna 2

DOA_angle_target=20;        % DOA of TARGET

sim_ant11=A_ant1*sin(2*pi*f*t);                         % Simulated antenna 1
sim_ant22=A_ant2*sin(2*pi*f*t+deg2rad(DOA_angle_target));   % Simulated antenna 2

sim_ant1=sim_ant11/max(sim_ant11);          %normalization
sim_ant2=sim_ant22/max(sim_ant22);

[smy,smx]=max(sim_ant1);

%%%%also take for same sample value of data

sim_sum_ant=sim_ant1(smx)+sim_ant2(smx);
sim_diff_ant=sim_ant1(smx)-sim_ant2(smx);
sim_ratio_ant=sim_diff_ant/sim_sum_ant;

%%%%%%%%%%%%% Polynomial fitting of ratio_ant of 6th order%%%%%%%%%%%%
% Error in DOA is obtained because of this
% if polynomial fitting is removed or more accurate results are obtained
% DOA of target can be found with greater accuracy
% Polynomial was obtained though curve fitting toolbox
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p1 =  1.047e-008;
p2 = -6.907e-007;
p3 =  2.383e-005;
p4 =  -0.0004914;
p5 =    0.008555;
p6 =    -0.09626;
p7 =      0.4215;
x=0;

for ii=1:2200
    x=x+0.01;
fitted_model_rat_ant(ii) = p1*x^6 + p2*x^5 + p3*x^4 + p4*x^3 + p5*x^2 + p6*x + p7;

end

theta_new=-scanagle_deg:.0545:scanagle_deg;    % Theta for ploting fitted model
figure
plot(theta_new(1:2200),fitted_model_rat_ant)

c1=ceil(sim_ratio_ant*7000);           % For comparison of sim_ratio ant and model
c2=ceil(fitted_model_rat_ant*7000);    % Different threshold Threshold can be chosen 
[r,c,v]= find(c2==c1);

detected_theta=(c.*0.0545)-60       %theta from curve fitting model

if(A_ant1>A_ant2)           % condition for checking which angle was correct
    correct_theta=detected_theta(1)
elseif(A_ant1<A_ant2)
    correct_theta=detected_theta(2)
else
     correct_theta=detected_theta
end

Correlation on DSP KIT TMS320C6713 Simulator

● December 15, 2010 ● Coded in C++ for the TI C67x
#include <stdio.h>
#include <math.h>

double pi=3.14159;

float x[50],y[50];
float r[100];

void correlation(float *x,float *y, int nx, int ny);

int main()
{   
   

int ii;
for(ii=0;ii<50  ;ii++)
 {
	 x[ii]=sin(2.0*pi*ii*3/50);
 	 y[ii]=sin(2.0*pi*ii*3/50);

 }

correlation(x,y,50,50);

printf("Complete.\n");

return 0;
}

void correlation(float *x,float *y, int nx, int ny)
{
	int n=50,delay=0,maxdelay=50;
	int i,j;
   double mx,my,sx,sy,sxy,denom;
  
   
   /* Calculate the mean of the two series x[], y[] */
   mx = 0;
   my = 0;   
	   for (i=0;i<n;i++) 
		   {
			  mx += x[i];
			  my += y[i];
		   }
   mx /= n;
   my /= n;

   /* Calculate the denominator */
   sx = 0;
   sy = 0;
	   for (i=0;i<n;i++) 
		   {
			  sx += (x[i] - mx) * (x[i] - mx);
			  sy += (y[i] - my) * (y[i] - my);
		   }
   denom = sqrt(sx*sy);

   /* Calculate the correlation series */
   for (delay=-maxdelay;delay<maxdelay;delay++) 
	   {
		  sxy = 0;
			  for (i=0;i<n;i++) 
				  {
					 j = i + delay;
	 if (j < 0 || j >= n)
						continue;
	 else
						sxy += (x[i] - mx) * (y[j] - my);
					 // Or should it be (?)
					/* if (j < 0 || j >= n)
						sxy += (x[i] - mx) * (-my);
					 else
						sxy += (x[i] - mx) * (y[j] - my);*/
					 
				  }
r[delay+maxdelay]= ( sxy / denom);

Linear and nonlinear phase low pass filter

● October 26, 2010●6 comments ● Coded in Matlab 
%% FIR & IIR filter
% Initializatoin

clc
clear all 
close all

%% FIR linear phase low pass filter

h=remez(60, [0 0.25 0.3 1], [1 1 0 0]);
fvtool(h,1)

%% IIR nonlinear phase low pass filter

[b a]=cheby2(20,35,0.3);
fvtool(b,a)

%% Filtering x(n)

x=zeros(1,200);
x(1:30)=1;

y1=filter(h,1,x);
y2=filter(b,a,x);

figure
plot(y1,'b','linewidth',2)
hold on
plot(y2,'r','linewidth',2)
xlabel('Samples')
ylabel('Amplitude')
title('Output of filters for input x(n)')
legend('FIR linear phase LPF output','IIR nonlinear phase LPF filter output')