Code Snippets Submitted by mdejesus
Discrete Fourier Transform (DFT)
function xdft=fdft(x)
N=length(x);
xdft(1,N)=0; % defining a 1xN sequence with all the elements equals to zero.
%*************THIS IS THE DFT IMPLEMENTATION**************
%For each value of k, the Summatory function is computed over the interval 0<n<N-1.
for k=0:N-1,
for n=0:N-1,
xdft(k+1)=xdft(k+1) + x(n+1)*exp(-j*2*pi*k*n/N);
end
end
%NOTE: You may skip the use of for-loops by defining k and n as two N-points sequences. This is:
% k=0:N-1; n=0:N-1;
%The DFT is now computed in matrix form (as a matrix multiplication of the N-points Fourier Matrix and the discrete-time sequence x).Goertzel Filterbank to the Implementation of a Nonuniform DFT
% VectorGoertzel Goertzel's Algorithm filter bank.
% Realization of the Goertzel's Algorithm to compute the Nonuniform DFT
% of a signal(a column vector named signalw) of length Nw with sampling
% frecuency fs at the desired frecuencies contained in vector f. The
% function returns the NDFT magnitude in a vector with the same length of f.
function xk=Gfilterbank(signalw,f,fs,Nw)
% Inititialization of the different variables
n=2;
signalw=signalw(1:Nw);
cost=cos(2*pi*f/fs);
sint=sin(2*pi*f/fs);
L=length(cost);
y1(1:L)=0;
y2(1:L)=0;
signalw=[0 0 signalw]; %Signal is delayed by two samples
% Goertzel Feedback Algorithm
while((n-2) < Nw)
n=n+1;
xnew(1:L)=signalw(n);
y=xnew+2*cost.*y1-y2;
y2=y1;
y1=y;
end
% Goertzel Forward Algorithm
rey=y1-y2.*cost;
imy=y2.*sint;
% Magnitude Calculation
xk=abs(rey+imy*j)';Beampattern of a Linear Array of Antennas (Array Processing)
% The function NDFT2 computes the sinc pattern of a linear array of
% sensors. The function receives three input: the interelement spacing d,
% the array weigths a, and the zero padding desired Npad.
function NDFT2(d, a, Npad)
k = -Npad/2 : Npad/2-1; % index
u = pi*k/Npad; % u is a vector with elements in the range of -pi/2 to pi/2
uk = sin(u);
n = 0:Npad-1; % n is an index
m = 1:Npad; % m is an index
% Wavenumber K = 2*pi/landa
% d = is a fraction or a multiple of lambda.
% v = K*d*uk(m).
v = 2*pi*d*uk(m)';
Dn(m,n+1) = exp(j*v*n);
puk = Dn * a'; % Computing the Beampattern
% Plot of the Beampatterns
figure(1); subplot(2,1,1); plot(uk,20*log10(abs(puk)));grid on; xlabel('sin(theta)'); ylabel('Magnitude (dB)')
axis([-1 1 -100 0])
subplot(2,1,2); plot(uk,puk);grid on; xlabel('sin(theta)'); ylabel('Magnitude')
axis([-1 1 -1 1])
warning off;
% Plot the beampattern as a polar graph
figure(2); polar(u',puk); hold on; polar(-u',puk);hold off
% Function pattern()
%
% The function pattern() computes and plots the beampattern of a
% Linear Array of sensors. This function has three inputs: the number of elements in
% the array, the pointing direction of the beam pattern (an angular
% direction in radians), and a constant spacing between the elements of
% the array (as fraction of the wavelenght(lambda)of the transmitted
% signal. The optimum spacing between elements of the array is 0.5.
% You must also select any of the windows presented in the menu.
% Windowing techniques are used to reduced the sidelobes of the pattern.
% The list of available windows is the following:
%
% @bartlett - Bartlett window.
% @barthannwin - Modified Bartlett-Hanning window.
% @blackman - Blackman window.
% @blackmanharris - Minimum 4-term Blackman-Harris window.
% @bohmanwin - Bohman window.
% @chebwin - Chebyshev window.
% @flattopwin - Flat Top window.
% @gausswin - Gaussian window.
% @hamming - Hamming window.
% @hann - Hann window.
% @kaiser - Kaiser window.
% @nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window.
% @parzenwin - Parzen (de la Valle-Poussin) window.
% @rectwin - Rectangular window.
% @tukeywin - Tukey window.
% @triang - Triangular window.
%
% For example: pattern(21, pi/4, 0.5);
% Plots the beampattern of an linear array with 21 elements equally spaced
% at a half of the wavelenght(lambda/2), and a pointing
% direction of 45 degrees. For uniform arrays use a rectangular
% window (rectwin).
function pattern(array_number, angular_direction, spacing)
close all
clc
N = array_number;
delta = angular_direction;
d = spacing;
Npad=1024;
n=0:N-1;
delta = 2*pi*d*sin(delta);
an=1/N*exp(j*n*delta);
for i=0:500000,
option = menu('Choose the desired Window', 'Bartlett', 'Barthannwin', 'Blackman', 'Blackmanharris', 'Bohmanwin', 'Chebwin', 'Flattopwin', 'Gausswin', 'Hamming', 'Hann', 'Kaiser', 'Nuttallwin', 'Parzenwin', 'Rectwin', 'Tukeywin', 'Triang', 'Exit');
switch option
case 1
close all
clear a;
a=an;
a = a.*bartlett(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 2
close all
clear a;
a=an;
a = a.*barthannwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 3
close all
clear a;
a=an;
a = a.*blackman(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 4
close all
clear a;
a=an;
a = a.*blackmanharris(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 5
close all
clear a;
a=an;
a = a.*bohmanwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 6
close all
clear a;
a=an;
a = a.*chebwin(N,40)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 7
close all
clear a;
a=an;
a = a.*flattopwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 8
close all
clear a;
a=an;
a = a.*gausswin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 9
close all
clear a;
a=an;
a = a.*hamming(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 10
close all
clear a;
a=an;
a = a.*hann(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 11
close all
clear a;
a=an;
a = a.*kaiser(N,1)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 12
close all
clear a;
a=an;
a = a.*nuttallwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 13
close all
clear a;
a=an;
a = a.*parzenwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 14
close all
clear a;
a=an;
a = a.*rectwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 15
close all
clear a;
a=an;
a = a.*tukeywin(N,0)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 16
close all
clear a;
a=an;
a = a.*triang(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 17
break;
end
endDiscrete Fourier Transform
#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
#define pi 3.14159265
// class complex
class complex
{
public:
//This class member function set the values of the class members real and imag
void setvals(double r, double i)
{
real=r;
imag=i;
}
//This function get the values (by reference) of the class members real and imag
void getvals(double &r, double &i)
{
r=real;
i=imag;
}
private:
/* The class members real and imag represents the real and imaginary parts
of a complex number, respectively.*/
double real;
double imag;
}; // end of class complex
//**********************************************************************************************
/* function dft receives two parameters x and N. x is the input sequence containing
the discrete-time samples and N is the number of DFT points. The complex type paramater
xk receives a sequence of N-points frequency samples (this is the DFT of the input sequence*/
void dft(double x[], int N, complex xk[])
{
double r,i;
for(int k=0;k<N;k++) //k is the index of the N-point frequency sequence
{
xk[k].setvals(0,0); //initialize the values of the counters to zero
for(int n=0;n<N;n++) //n is the index of the discrete-time sequence
{
xk[k].getvals(r,i); // get the values of real and imag
/*This is the computation of the real and imaginary parts of the DFT.
The class member function setvals is used to update the value of the accumulators
contaning the real and imaginary parts of the DFT samples*/
xk[k].setvals(r + x[n]*cos(2*pi*k/N*n),i + x[n]*sin(2*pi*k/N*n));
} //end inner for-loop
} //end outer for-loop
}//end dft function
//************************************************************************************************
//Declaration of the main() function
int main()
{
int size; // The number of DFT samples
double r,i; // r = real part of the DFT; i = imaginary part of the DFT
cout<<"Enter the desired number of DFT samples: ";
cin>>size;
double *x = new double[size]; //dynamic array x represent the discrete-time samples
complex *xk = new complex[size]; // dynamic array xk represent the DFT samples
cout<<"\nEnter the sequence x[n]: ";
for(int h=0;h<size;h++)
{
cin>>x[h];
}
//Computation of the DFT of x
dft(x, size, xk);
system("CLS"); //clear screen
cout<<"Discrete Fourier Transform: "<<endl<<endl;
for(int h=0;h<size;h++)
{
cout.precision(3); // display three (3) decimal places
xk[h].getvals(r,i); //get the DFT samples values
cout<<"xk["<<h<<"] = "<<fixed<< r <<" + "<<fixed<< i <<"j"<<endl; // print the DFT samples
}
cout<<endl;
return 0;
}// end main functionEcho Filter
% function echo_filter:
% 1) computes and returns the impulse response h of the LTI system (echo filter)
% 2) Filter the input signal. The output signal is yecho.
%
% INPUTS:
% signal = input signal
% times = vector containing the time delays (in seconds) of the repetitions(echoes)
% attenuations = vector containing the respective attenuations of the echoes
% fs = sampling frequency of the input signal.
%
% OUTPUTS:
% yecho = output signal (passed through the echo filter)
% h = impulse response of the echo filter.
function [yecho, h] = echo_filter(signal, times, attenuations, fs)
h(1)=1; % This is the first coefficient of the echo filter
% Computing the impulse response h
for i=1:length(times),
samples(i) = times(i)*fs; % Calculating the sample-times (N = t*fs)
h(floor(samples(i))) = attenuations(i); % Impulse response coeficients
end
% #########################################################################
% You may use the following implementation instead of the illustrated
% above:
% samples = times*fs;
% h(floor(samples)) = attenuations;
% In this case, the implementation is done without a for-loop.
% #########################################################################
yecho = filter(h,1,signal(:,1)); % filtering of the input signal results in a signal with echoes
% You may test the filter using the Matlab signal "splat" (write load splat
% in the Matlab command window). Use function sound() in order to
% listening both, the input and the output, signals.Goertzel Filterbank to the Implementation of a Nonuniform DFT
% VectorGoertzel Goertzel's Algorithm filter bank.
% Realization of the Goertzel's Algorithm to compute the Nonuniform DFT
% of a signal(a column vector named signalw) of length Nw with sampling
% frecuency fs at the desired frecuencies contained in vector f. The
% function returns the NDFT magnitude in a vector with the same length of f.
function xk=Gfilterbank(signalw,f,fs,Nw)
% Inititialization of the different variables
n=2;
signalw=signalw(1:Nw);
cost=cos(2*pi*f/fs);
sint=sin(2*pi*f/fs);
L=length(cost);
y1(1:L)=0;
y2(1:L)=0;
signalw=[0 0 signalw]; %Signal is delayed by two samples
% Goertzel Feedback Algorithm
while((n-2) < Nw)
n=n+1;
xnew(1:L)=signalw(n);
y=xnew+2*cost.*y1-y2;
y2=y1;
y1=y;
end
% Goertzel Forward Algorithm
rey=y1-y2.*cost;
imy=y2.*sint;
% Magnitude Calculation
xk=abs(rey+imy*j)';Discrete Fourier Transform
#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
#define pi 3.14159265
// class complex
class complex
{
public:
//This class member function set the values of the class members real and imag
void setvals(double r, double i)
{
real=r;
imag=i;
}
//This function get the values (by reference) of the class members real and imag
void getvals(double &r, double &i)
{
r=real;
i=imag;
}
private:
/* The class members real and imag represents the real and imaginary parts
of a complex number, respectively.*/
double real;
double imag;
}; // end of class complex
//**********************************************************************************************
/* function dft receives two parameters x and N. x is the input sequence containing
the discrete-time samples and N is the number of DFT points. The complex type paramater
xk receives a sequence of N-points frequency samples (this is the DFT of the input sequence*/
void dft(double x[], int N, complex xk[])
{
double r,i;
for(int k=0;k<N;k++) //k is the index of the N-point frequency sequence
{
xk[k].setvals(0,0); //initialize the values of the counters to zero
for(int n=0;n<N;n++) //n is the index of the discrete-time sequence
{
xk[k].getvals(r,i); // get the values of real and imag
/*This is the computation of the real and imaginary parts of the DFT.
The class member function setvals is used to update the value of the accumulators
contaning the real and imaginary parts of the DFT samples*/
xk[k].setvals(r + x[n]*cos(2*pi*k/N*n),i + x[n]*sin(2*pi*k/N*n));
} //end inner for-loop
} //end outer for-loop
}//end dft function
//************************************************************************************************
//Declaration of the main() function
int main()
{
int size; // The number of DFT samples
double r,i; // r = real part of the DFT; i = imaginary part of the DFT
cout<<"Enter the desired number of DFT samples: ";
cin>>size;
double *x = new double[size]; //dynamic array x represent the discrete-time samples
complex *xk = new complex[size]; // dynamic array xk represent the DFT samples
cout<<"\nEnter the sequence x[n]: ";
for(int h=0;h<size;h++)
{
cin>>x[h];
}
//Computation of the DFT of x
dft(x, size, xk);
system("CLS"); //clear screen
cout<<"Discrete Fourier Transform: "<<endl<<endl;
for(int h=0;h<size;h++)
{
cout.precision(3); // display three (3) decimal places
xk[h].getvals(r,i); //get the DFT samples values
cout<<"xk["<<h<<"] = "<<fixed<< r <<" + "<<fixed<< i <<"j"<<endl; // print the DFT samples
}
cout<<endl;
return 0;
}// end main functionEcho Filter
% function echo_filter:
% 1) computes and returns the impulse response h of the LTI system (echo filter)
% 2) Filter the input signal. The output signal is yecho.
%
% INPUTS:
% signal = input signal
% times = vector containing the time delays (in seconds) of the repetitions(echoes)
% attenuations = vector containing the respective attenuations of the echoes
% fs = sampling frequency of the input signal.
%
% OUTPUTS:
% yecho = output signal (passed through the echo filter)
% h = impulse response of the echo filter.
function [yecho, h] = echo_filter(signal, times, attenuations, fs)
h(1)=1; % This is the first coefficient of the echo filter
% Computing the impulse response h
for i=1:length(times),
samples(i) = times(i)*fs; % Calculating the sample-times (N = t*fs)
h(floor(samples(i))) = attenuations(i); % Impulse response coeficients
end
% #########################################################################
% You may use the following implementation instead of the illustrated
% above:
% samples = times*fs;
% h(floor(samples)) = attenuations;
% In this case, the implementation is done without a for-loop.
% #########################################################################
yecho = filter(h,1,signal(:,1)); % filtering of the input signal results in a signal with echoes
% You may test the filter using the Matlab signal "splat" (write load splat
% in the Matlab command window). Use function sound() in order to
% listening both, the input and the output, signals.Beampattern of a Linear Array of Antennas (Array Processing)
% The function NDFT2 computes the sinc pattern of a linear array of
% sensors. The function receives three input: the interelement spacing d,
% the array weigths a, and the zero padding desired Npad.
function NDFT2(d, a, Npad)
k = -Npad/2 : Npad/2-1; % index
u = pi*k/Npad; % u is a vector with elements in the range of -pi/2 to pi/2
uk = sin(u);
n = 0:Npad-1; % n is an index
m = 1:Npad; % m is an index
% Wavenumber K = 2*pi/landa
% d = is a fraction or a multiple of lambda.
% v = K*d*uk(m).
v = 2*pi*d*uk(m)';
Dn(m,n+1) = exp(j*v*n);
puk = Dn * a'; % Computing the Beampattern
% Plot of the Beampatterns
figure(1); subplot(2,1,1); plot(uk,20*log10(abs(puk)));grid on; xlabel('sin(theta)'); ylabel('Magnitude (dB)')
axis([-1 1 -100 0])
subplot(2,1,2); plot(uk,puk);grid on; xlabel('sin(theta)'); ylabel('Magnitude')
axis([-1 1 -1 1])
warning off;
% Plot the beampattern as a polar graph
figure(2); polar(u',puk); hold on; polar(-u',puk);hold off
% Function pattern()
%
% The function pattern() computes and plots the beampattern of a
% Linear Array of sensors. This function has three inputs: the number of elements in
% the array, the pointing direction of the beam pattern (an angular
% direction in radians), and a constant spacing between the elements of
% the array (as fraction of the wavelenght(lambda)of the transmitted
% signal. The optimum spacing between elements of the array is 0.5.
% You must also select any of the windows presented in the menu.
% Windowing techniques are used to reduced the sidelobes of the pattern.
% The list of available windows is the following:
%
% @bartlett - Bartlett window.
% @barthannwin - Modified Bartlett-Hanning window.
% @blackman - Blackman window.
% @blackmanharris - Minimum 4-term Blackman-Harris window.
% @bohmanwin - Bohman window.
% @chebwin - Chebyshev window.
% @flattopwin - Flat Top window.
% @gausswin - Gaussian window.
% @hamming - Hamming window.
% @hann - Hann window.
% @kaiser - Kaiser window.
% @nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window.
% @parzenwin - Parzen (de la Valle-Poussin) window.
% @rectwin - Rectangular window.
% @tukeywin - Tukey window.
% @triang - Triangular window.
%
% For example: pattern(21, pi/4, 0.5);
% Plots the beampattern of an linear array with 21 elements equally spaced
% at a half of the wavelenght(lambda/2), and a pointing
% direction of 45 degrees. For uniform arrays use a rectangular
% window (rectwin).
function pattern(array_number, angular_direction, spacing)
close all
clc
N = array_number;
delta = angular_direction;
d = spacing;
Npad=1024;
n=0:N-1;
delta = 2*pi*d*sin(delta);
an=1/N*exp(j*n*delta);
for i=0:500000,
option = menu('Choose the desired Window', 'Bartlett', 'Barthannwin', 'Blackman', 'Blackmanharris', 'Bohmanwin', 'Chebwin', 'Flattopwin', 'Gausswin', 'Hamming', 'Hann', 'Kaiser', 'Nuttallwin', 'Parzenwin', 'Rectwin', 'Tukeywin', 'Triang', 'Exit');
switch option
case 1
close all
clear a;
a=an;
a = a.*bartlett(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 2
close all
clear a;
a=an;
a = a.*barthannwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 3
close all
clear a;
a=an;
a = a.*blackman(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 4
close all
clear a;
a=an;
a = a.*blackmanharris(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 5
close all
clear a;
a=an;
a = a.*bohmanwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 6
close all
clear a;
a=an;
a = a.*chebwin(N,40)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 7
close all
clear a;
a=an;
a = a.*flattopwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 8
close all
clear a;
a=an;
a = a.*gausswin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 9
close all
clear a;
a=an;
a = a.*hamming(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 10
close all
clear a;
a=an;
a = a.*hann(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 11
close all
clear a;
a=an;
a = a.*kaiser(N,1)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 12
close all
clear a;
a=an;
a = a.*nuttallwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 13
close all
clear a;
a=an;
a = a.*parzenwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 14
close all
clear a;
a=an;
a = a.*rectwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 15
close all
clear a;
a=an;
a = a.*tukeywin(N,0)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 16
close all
clear a;
a=an;
a = a.*triang(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 17
break;
end
endDiscrete Fourier Transform (DFT)
function xdft=fdft(x)
N=length(x);
xdft(1,N)=0; % defining a 1xN sequence with all the elements equals to zero.
%*************THIS IS THE DFT IMPLEMENTATION**************
%For each value of k, the Summatory function is computed over the interval 0<n<N-1.
for k=0:N-1,
for n=0:N-1,
xdft(k+1)=xdft(k+1) + x(n+1)*exp(-j*2*pi*k*n/N);
end
end
%NOTE: You may skip the use of for-loops by defining k and n as two N-points sequences. This is:
% k=0:N-1; n=0:N-1;
%The DFT is now computed in matrix form (as a matrix multiplication of the N-points Fourier Matrix and the discrete-time sequence x).
