DSPRelated.com
Code
   Code Snippets    Submitted by mdejesus

Goertzel Filterbank to the Implementation of a Nonuniform DFT

Miguel De Jesus Rosario ● December 14, 2010●1 comment ● Coded in Matlab 
% 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

Miguel De Jesus Rosario ● December 14, 2010●3 comments ● Coded in C++ 
#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 function

Echo Filter

Miguel De Jesus Rosario ● December 1, 2010●3 comments ● Coded in Matlab 
% 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)

Miguel De Jesus Rosario ● November 30, 2010●1 comment ● Coded in Matlab 
% 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
    
end

Discrete Fourier Transform (DFT)

Miguel De Jesus Rosario ● October 27, 2010●11 comments ● Coded in Matlab 
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).