DSPRelated.com
   Code Snippets

Peak/Notch Filter Design

Jeff T ● May 21, 2011●3 comments ● Coded in Matlab 
function [b, a]  = peaking(G, fc, Q, fs)

% Derive coefficients for a peaking filter with a given amplitude and
% bandwidth.  All coefficients are calculated as described in Zolzer's
% DAFX book (p. 50 - 55).  This algorithm assumes a constant Q-term
% is used through the equation.
%
% Usage:     [B,A] = peaking(G, Fc, Q, Fs);
%
%            G is the logrithmic gain (in dB)
%            FC is the center frequency
%            Q is Q-term equating to (Fb / Fc)
%            Fs is the sampling rate
%
% Author:    sparafucile17 08/22/05
%

K = tan((pi * fc)/fs);
V0 = 10^(G/20);

%Invert gain if a cut
if(V0 < 1)
    V0 = 1/V0;
end

%%%%%%%%%%%%%%
%   BOOST
%%%%%%%%%%%%%%
if( G > 0 )
   
    b0 = (1 + ((V0/Q)*K) + K^2) / (1 + ((1/Q)*K) + K^2);
    b1 =        (2 * (K^2 - 1)) / (1 + ((1/Q)*K) + K^2);
    b2 = (1 - ((V0/Q)*K) + K^2) / (1 + ((1/Q)*K) + K^2);
    a1 = b1;
    a2 =  (1 - ((1/Q)*K) + K^2) / (1 + ((1/Q)*K) + K^2);
  
%%%%%%%%%%%%%%
%    CUT
%%%%%%%%%%%%%%
else
    
    b0 = (1 + ((1/Q)*K) + K^2) / (1 + ((V0/Q)*K) + K^2);
    b1 =       (2 * (K^2 - 1)) / (1 + ((V0/Q)*K) + K^2);
    b2 = (1 - ((1/Q)*K) + K^2) / (1 + ((V0/Q)*K) + K^2);
    a1 = b1;
    a2 = (1 - ((V0/Q)*K) + K^2) / (1 + ((V0/Q)*K) + K^2);
    
end

%return values
a = [  1, a1, a2];
b = [ b0, b1, b2];

Interpolation from Nonuniform Samples

Jesus Selva ● May 18, 2011●1 comment ● Coded in Matlab 
function Demo

%Author: J. Selva.
%Date: May 2011.

%The interpolation method in this snippet has been published in 
%
% [1] J. Selva, "Functionally weighted Lagrange interpolation of bandlimited
% signals from nonuniform samples," IEEE Transactions on Signal Processing, 
% vol. 57, no. 1, pp. 168-181, Jan 2009.
%
% Also, the barycentric implementation of the interpolator appears in 
%
% [2] J. Selva, "Design of barycentric interpolator for uniform and nonuniform 
% sampling grids", IEEE Trans. on Signal Processing, vol. 58, n. 3, 
% pp. 1618-1627, March 2010.

T = 1;   %Sampling period
B = 0.7; %Two-sided signal bandwidth. It must be B*T < 1.
P = 12;  %Truncation index. The number of samples taken is 2*P+1. The error 
         %decreases EXPONENTIALLY with P as exp(-pi*(1-B*T)*P). 
 	%(Try different values of P like 3, 7, 20, and 30.)

int = [-30,30]*T; %Plot x limits.

JitMax = 0.4; %Maximum jitter. This must be between 0 and 0.5*T. The samples
              %are taken at instants p*T+d[p], -P<=p<=P and 
	      %-JitMax <= d[p] <= JitMax.  

while 1
  
  st = TestSignal(int,ceil(60/B),B); %This defines a test signal formed 
              %by a random sum of sincs of bandwidth B. 
  
  d = 2*(rand(1,2*P+1)-0.5)*JitMax; %Vector of random jitter.

  disp('Sampling instants (t_k/T):')
  disp((-P:P)*T+d)
  
  sg = TestSignal(st,(-P:P)*T+d); %Sample signal at jittered instants

  t = int(1):0.1*T:int(2); %Time grid for figure.

  sRef = TestSignal(st,t); %Exact signal samples in grid t.
  sInt = NonuniformInterp(t,sg,d,B,T); %Signal samples interpolated from 
                           %d (jitter values) and sg (values at jittered 
			   %instants).

  subplot(2,1,1)
  plot(t,sRef,t,sInt) %Plot exact and interpolated signals

  xlabel('Normalized time (t/T)')
  title('Signal (blue) and interpolator (green)')
  grid on

  subplot(2,1,2)
  plot(t,20*log10(abs(sRef-sInt))) %Plot interpolation error in dB.
  xlabel('Normalized time (t/T)')
  ylabel('Interpolation error (dB)')
  grid on
  
  R = input('Press Enter for a new trial, q to quit: ','s');

  if any(strcmp(R,{'q','Q'}))
    break
  end

end

function v = NonuniformInterp(u,sg,d,B,T)

%Performs interpolation from nonuniform samples.
% T   is the period of the reference grid (instants (-P:P)*T).
% d   is the vector of jitter values relative to the reference grid. It must have an odd ...
%     number of elements. 
% B:  signal's two-sided bandwidth. It must be B*T<1. The error decreases exponentially as
%     exp(-pi*(1-B*T)P).
% u:  Vector of instants where the signal is to be interpolated.
% sg: Sample values at instants (-P:P)*T+d.

if mod(length(d),2) ~= 1
  error('The number of samples must be odd')
end

P = (length(d)-1)/2;

tg = (-P:P)*T+d(:).';

w = baryWeights(tg,T,B,P);

v =  DerBarycentricInterp3tVecV1(w,sg,tg,u,0);

function v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function w = baryWeights(vt,T,B,P)

%Barycentric weights. See Eq. (23) and sequel in [2].

g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;

N = length(vt);
LD = ones(1,N);

for k = 1:N
  LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end

w = gam ./ LD;
w = w / max(abs(w));

function out = DerBarycentricInterp3tVecV1(alpha,s,t,tau,nd)

%Vectorized implementation of barycentric interpolator in [2, Sec. IV].

s = s(:).';
t = t(:).';
sztau = size(tau);
tau = tau(:);

Nt = numel(t);
Ntau = numel(tau);

vD = zeros(Ntau,1);
LF = [ones(Ntau,1),zeros(Ntau,Nt-1)];
out = zeros(Ntau,nd+1);

for k = 1:Nt-1
  vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
  LF(:,k+1) = LF(:,k) .* (tau-t(k));
end

vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);

z = s(ones(Ntau,1),:);

for kd = 0:nd

  pr = z(:,end);
  z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);

  for k = 1:Nt-1
    cN = cN .* (tau-t(k))+z1(:,k)*alpha(k) .* LF(:,k);
  end
  cN = cN ./ vD;
  
  ztau = z(:,end)+(tau-t(end)) .* cN;
  out(:,kd+1) = ztau;
  
  if kd < nd
    z = [ (z(:,1:end-1)-ztau)./(t(ones(Nt,1),1:end-1)-tau(:,ones(1,Nt))) , ...
	  cN ] * (kd+1);
  end
  
end

out = reshape(out,sztau);

function varargout = TestSignal(varargin)

%Test signal formed by a random sum of sincs. 

if nargin == 3
  [int,Nsinc,B] = deal(varargin{:});
  st = [];
  st.B = B;
  st.Ampl = (rand(1,Nsinc)-0.5)*2;
  st.Del = int(1)+(int(2)-int(1))*rand(1,Nsinc);
  varargout = {st};
else
  [st,t] = deal(varargin{:});
  v = zeros(size(t));
  for k = 1:numel(t)
    v(k) = st.Ampl * sinc(st.B*(t(k)-st.Del)).';
  end
  varargout = {v};
end

T1 F281x Mcbsp Transmitter Configuration

Emmanuel ● May 13, 2011 ● Coded in C for the TI F28x
/**********************************************************************/
#include "DSP281x_Device.h"
/**********************************************************************
Parameter					Condition		Register		Value

Delay between data frames:	1 bit			XCR2.XDATDLY	0b01
Word Length:				16 bits			XCR1.XDLEN1		0b010
Words per frame:			1 word			XCR1.XFRLEN1	0
sample rate:				150 kHz			SRGR1.CLKGDV	0xC7
															(LSPCLK=30 MHz)

************************************************************************/
void InitMcbsp(void)
{

/*** Mcbsp	Reset	*/

	McbspaRegs.SPCR2.bit.XRST= 0;
	McbspaRegs.SPCR2.bit.GRST= 0;
	McbspaRegs.SPCR2.bit.FRST= 0;
	
/*** SPCR2 Register Configuration	*/

	McbspaRegs.SPCR2.all = 0x0000;
	
// bit 15-10 	0's:	Reserved
// bit 9 		0:		FREE Free running mode
// bit 8 		0:		SOFT Soft bit
// bit 7 		0:		FRST Frame sync generator reset
// bit 6 		0:		GRST Sample rate generator reset
// bit 5-4 		00:		XINTM Transmit interrupt mode
// bit 3 		0:		XSYNCERR Transmit synchronization error
// bit 2 		0:		XEMPTY Transmit shift register (XSR) empty
// bit 1 		0:		XRDY Transmitter ready
// bit 0 		0:		XRST Transmitter reset
	
/*** SPCR1 Register Configuration 		*/	

	McbspaRegs.SPCR1.all = 0x0000;

// bit 15 		0:		DLB Digital Loop back mode
// bit 14-13 	0:		RJUST Receive sign extension and justification mode
// bit 12-11 	00:		CLKSTP Clock Stop Mode
// bit 10-8 	0's:	Reserved
// bit 7 		0:		DXENA DX enabler. 
// bit 6 		0:		ABIS ABIS mode
// bit 5-4 		00:		RINTM Receive Interrupt Mode
// bit 3 		0		RSYNCERR Receive Synchronization error
// bit 2 		0		RFULL Receive shift register full
// bit 1 		0		RRDY Receiver ready
// bit 0 		0		RRST Receiver reset

/*** XCR2 Register Configuration 	*/	

	McbspaRegs.XCR2.all = 0x0001;

// bit 15 		0:		XPHASE Transmit phase
// bit 14-8 	0's;	XRFRLEN2 Transmit frame length 2
// bit 7-5 		000:	XWDLEN2 Transmit word length 2
// bit 4-3 		00:		XCOMPAND Transmit companding mode
// bit 2 		0:		XFIG Transmite frame ignore
// bit 1-0 		01:		XDATDLY Data Delay (Delay in bits between two frames)
	
/*** XCR1 Register Configuration 	*/	

	McbspaRegs.XCR1.all = 0x0040;
	
// bit 15 		0		Reserved
// bit 14-8 	0's:	XFRLEN1 Frame Length 1 (Words per frame)
// bit 7-5 		010:	XWDLEN1 Word Length 1  (bits per word)
// bit 4 		0:		Reserved
// bit 3-0 		0's:	Reserved

/*** SRGR2 Register Configuration 	*/

	McbspaRegs.SRGR2.all = 0x2011;

// bit 15 		0;		GSYNC Sample rate generator clock synchronization
// bit 14 		0:		Reserved
// bit 13 		1:		CLKSM Mcbsp sample rate generator clock mode: internal clock
// bit 12 		0:		FSGM Transmit frame synchronization mode
// bit 11-0 	x011:	FPER Frame Period (number of CLKG cycles between sync signals): 18 cycles

/*** SRGR1 Register Configuration 	*/	

	McbspaRegs.SRGR1.all = 0x00E8;

// bit 15-8 	x00:	FWID Frame Width
// bit 7-0 		xC7		CLKGDV Sample rate generator clock divider

/*** MCR2 Register Configuration 	*/	

	McbspaRegs.MCR2.all = 0x0200;
	
// bit 15-10 	0's:	Reserved
// bit 9 		1:		XMCME Enhanced transmit multichannel selection enable
// bit 8-7 		00:		XPBBLK Receive/transmit partition B block
// bit 6-5 		00:		XPABLK Receive/transmit partition A block
// bit 4-2 		000:	XCBLK Receive/transmit current block
// bit 1-0 		00:		XMCM Transmit multichannel selection enable

/*** MCR1 Register Configuration  */	

	McbspaRegs.MCR1.all = 0x0201;
	
// bit 15-10 	0's:	Reserved
// bit 9 		1		RMCME Enhanced receive multichannel selection enable
// bit 8-7 		00		RPBBLK Receive/transmit partition B block
// bit 6-5 		00		RPABLK Receive/transmit partition A block
// bit 4-2 		000		RCBLK Receive transmit current block
// bit 1 		0		Reserved
// bit 0 		1		RMCM Receive multichannel selection enable

/*** PCR Register Configuration */	

	McbspaRegs.PCR.all = 0x0A00;

// bit 15-12 	0's:	Reserved
// bit 11 		1:		FSXM Transmit frame synchronization mode
// bit 10 		0:		FSRM Receive frame synchronization mode
// bit 9 		1:		CLKXM Transmit frame synchronization mode (internal, FSR as output)
// bit 8 		0:		CLKRM Receiver clock mode
// bit 7 		0:		SCLKME Sample clock mode selection bit
// bit 6 		0:		CLKS_STAT Reserved 
// bit 5 		0:		DX_STAT DX pin status
// bit 4 		0:		DR_STAT DR pin status
// bit 3 		0:		FSXP Transmit frame synchronization polarity
// bit 2 		0:		FSRP Receive frame synchronization polarity
// bit 1 		0:		CLKXP Tramsmit clock polarity
// bit 0 		0:		CLKRP Receive clock polarity
	

	
/*** Mcbsp Start-up		*/

	McbspaRegs.SPCR2.bit.XRST= 1;
	McbspaRegs.SPCR2.bit.GRST= 1;
	McbspaRegs.SPCR2.bit.FRST= 1;
	

/*** Interruption Configuration  */

	McbspaRegs.MFFINT.bit.XINT= 1;

	PieCtrlRegs.PIEIER6.bit.INTx6= 1;	
    IER |= 0x0020;   
    
    

}

/**End of function*/

Flanger Audio Effect

Gabriel Rivas ● April 22, 2011 ● Coded in C 
/*******************FLANGER.C******************************/

#include "Delay.h"
#include "Flanger.h"

static short samp_freq;
static double var_delay;
static short counter;
static short counter_limit;
static short control;
static short max_delay;
static short min_delay;
static double mix_vol;
static double delay_step;

/*
This is the initialization function, basically
it passes the initialization parameters to the delay block
and initializes the flanger control variables.
*/
void Flanger_init(short effect_rate,short sampling,short maxd,short mind,double fwv,double stepd,double fbv) {
	Delay_Init(2,fbv,fwv,1);
 
	samp_freq = sampling;
	counter = effect_rate;
	control = 1;
	var_delay = mind;

	//User Parameters
	counter_limit = effect_rate;
	max_delay =  maxd;
	min_delay = mind;
	mix_vol = 1;
	delay_step = stepd;
}

/*This is the flanging process task
that uses the delay task inside*/
double Flanger_process(double xin) {
	double yout;

	yout = Delay_task(xin); 	
	return yout;
}

/*
This sweep function creates a slow frequency
ramp that will go up and down changing the
delay value at the same time. The counter
variable is a counter of amount of samples that 
the function waits before it can change the delay.
*/
void Flanger_sweep(void) {
	if (!--counter) {			
	    var_delay+=control*delay_step;
  
	    if (var_delay > max_delay) {
	        control = -1;
            } 

            if (var_delay < min_delay) {
	        control = 1;
            }

	    Delay_set_delay(var_delay);
	    counter = counter_limit;
	}
}

/*****************FLANGER.H********************************/
#ifndef __FLANGER_H__
#define __FLANGER_H__

extern void Flanger_init(short effect_rate,short sampling,short maxd,short mind,double fwv,double stepd,double fbv);
extern double Flanger_process(double xin);
extern void Flanger_sweep(void);

#endif

/*****USAGE EXAMPLE****************************************/
#include "Flanger.h"

void main(void) {
    double xin;
    double yout;
    Flanger_init(500,16000,70,2,0.3,1,0.3);

    while(1) {
        if (new_sample_flag()) {
            /*When there's new sample at your ADC or CODEC input*/
            /*Read the sample*/
            xin = read_sample();
            /*Apply the Flanger_process function to the sample*/
            yout = Flanger_process(0.7*xin);

            /*Send the output value to your DAC or codec output*/
            write_output(yout);
            /*Makes the delay vary*/
            Flanger_sweep();
        }                              
    }
}

TI F281x PWM Configuration

Emmanuel ● April 18, 2011 ● Coded in C for the TI F28x
/**********************************************************************/
#include "DSP281x_Device.h"
/***********************************************************************/

void InitPwm(void)
{

asm( "EALLOW");	//Enable write operations in protected registers

/*** Independent Compare Output Enable	*/

	EvaRegs.EXTCONA.bit.INDCOE= 1;

/********************************************************************
*	Write values in the two following registers until all the other
*	registers have been written, as their value depends on the other
*	registers. 							*
********************************************************************/

/*** Timer 1 Period Regiser	*/

	EvaRegs.T1PR = 0xB71A;			//use this value for a frequency of 200 Hz assuming
								    //sysclkout=150 MHz and its pre-scaler=16
								    
/*** CMPR1 Register Configuration */

	EvaRegs.CMPR1= 0x493D;		//to get 60% of duty cycle this value equals the 1-0.6= 40% period
		
/********************************************************************/

	EvaRegs.GPTCONA.all= 0x0050;
	
// bit 15 		0		Reserved 
// bit 14 		0		T2STAT GP timer 2 Status. 
// bit 13 		0		T1STAT GP timer 1 Status. 
// bit 12 		0		T2CTRIPE T2CTRIP Enable. 
// bit 11 		0		T1CTRIPE T1CTRIP Enable. 
// bit 10-9 	00		T2TOADC Start ADC with timer 2 event
// bit 8-7 		00		T1TOADC Start ADC with timer 1 event
// bit 6 		1		TCMPOE Timer compare output enable. 
// bit 5		0		T2CMPOE Timer 2 compare output enable. 
// bit 4 		1		T1CMPOE Timer 1 Compare Output Enable. 
// bit 3-2 		00		T2PIN Polarity of GP timer 2 compare output
// bit 1-0 		00		T1PIN Polarity of GP timer 1 compare output
	

/*** COMCONA Register Configuration*/

	EvaRegs.COMCONA.all = 0x8020;
	
// bit 15 		1:		CENABLE: Compare Enable
// bit 14-13 	00:		CLD1, CLD0: CMPRx Reload Condition
// bit 12 		0:		SVENABLE: PWM Vector Enable
// bit 11-10 	00:		ACTRLD1,ACTRLD0: Action Control Regiser Reload Condition
// bit 9 		0:		FCMPOE: Full Compare Output Enable
// bit 8 		0:		PDPINTA: PDPINTA Status
// bit 7 		0:		FCMP3OE: Compare 3 Output Enable
// bit 6 		0:		FCMP2OE: Compare 2 Output Enable
// bit 5 		1:		FCMP1OE: Compare 1 Output Enable
// bit 4-3 		00:		Reserved
// bit 2 		0:		C3TRIPE: C3TRIP Enable 
// bit 1 		0:		C2TRIPE: C2TRIP Enable 
// bit 0 		0:		C1TRIPE: C1TRIP Enable 

/** ACTRA Register Configuration*/

	EvaRegs.ACTRA.all = 0x0002;
	
// bit 15 		0:		SVRDIR: PWM Vector rotate direction 
// bit 14-12 	000:	D2-D: PWM Vector basic bits
// bit 11-10 	00:		CMP6ACT1-0 Compare Action in pin 6, CMP6.
// bit 9-8 		00:		CMP5ACT1-0 Compare Action in pin 5, CMP5.
// bit 7-6 		00:		CMP4ACT1-0 Compare Action in pin 4, CMP4.
// bit 5-4 		00:		CMP3ACT1-0 Compare Action in pin 3, CMP3
// bit 3-2 		00:		CMP2ACT1-0 Compare Action in pin 2, CMP2
// bit 1-0 		10		CMP1ACT1-0 Compare Action in el pin 1, CMP1

/*** T1CON Register Configuration*/

	EvaRegs.T1CON.all = 0x9442;
	
// bit 15:14 	10:		Emulation Control bits
// bit 13 		0:		Reserved
// bit 12-11 	10:		TMODE1-TMODE0 Counter Mode Selection
// bit 10-8 	100:	TPS2-TPS0 Pre-scaler: 100= HSPCLK/16
// bit 7 		0:		T2SWT1 Timer 2 trigger by timer 1
// bit 6 		1:		TENABLE Timer Enable
// bit 5-4 		00:		TCLKS(1,0) Clock Source Selection
// bit 3-2 		00:		TCLD(1,0) Timer Compare Register Reload Condition
// bit 1 		1:		TECMPR Compare Operations Enable
// bit 0		0:		SELT1PR: Period Register Selection
						

asm( "EDIS");	//Protected Registers Write Disabled

}

/**End of function*/

ML frequency estimation using the FFT

Jesus Selva ● April 18, 2011 ● Coded in Matlab 
function Demo(L)

%Author: J. Selva. 
%Date: April 2011.

%This demo constructs trials of an undamped exponential plus noise, and then computes the
%maximum likelihood estimate of its frequency, amplitude, and phase. Only one FFT of size
%2^ceil(log2(st.M)+1) is computed, plus a few operations. The complexity is that given by
%the FFT. The algoritm is explained in 
%
%         http://arxiv.org/abs/1104.3069
%

format long g
format compact

M = 1024; %Number of samples.
gdB = 20; %SNR at matched filter output.
ARef = 1.3; %Amplitude.
fRef = 0.372; %Frequency in [-0.5,0.5[.
PhiRef = 0.43; %Phase (rad).

if nargin > 0
  for k = 1:2:length(L)
    eval([L{k} '=L{k+1};']);
  end
end

%Define a random vector with 1024 elements. 

SNRdB =  gdB - 10*log10(M); %Signal-to-noise ratio (dB).

%Complex sinusoidal amplitude, frequency, and phase.

NVar = ARef^2*10.^(-SNRdB/10); %Noise variance. 

m1 = -ceil(M/2);

while 1
  
  disp(' ')
  disp(['M = ', num2str(M), ', gdB = ', num2str(gdB), ' dB, ARef = ', num2str(ARef), ...
	', fRef = ', num2str(fRef), ', PhiRef = ', num2str(PhiRef), ' rad'])
  disp(' ')
  
  %Set up a trial.
  a = ARef*exp(1i*2*pi*fRef*(m1:m1+M-1).'+1i*PhiRef)+(randn(M,1)+1i*randn(M,1))*sqrt(NVar/2);

  %Call the nufft function. The output is a struct with some variables for later use.

  st = nufft(a);

  %Compute the DFT peak. 
  [f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st);

  %The output is the following, 
  %
  %   f1, A1, and Phi1 are the ML estimates of frequency, amplitude and phase, respectively. 
  %
  %   L1 is the ML cost function value at f1 and its first two  derivatives. L1(1) should be
  %     large   since it is   the square value of   the peak amplitude.    L1(2) is the cost
  %     function   derivative,  and  it  should  be   small  since  we  are  at   a critical
  %     point. Finally  L1(2) should be a  negative number with large  amplitude since f1 is
  %     the abcissa of a maximum (second derivative).
  %
  %   Deg: If Deg==1 then the Newton iteration was cut too many times. 
  %   nIt: Number of Newton iterations. 

  disp('Exact and estimated frequency, amplitude, and phase:')
  disp([fRef,ARef,PhiRef;f1,A1,Phi1])
  disp(['Number of iterations: ' num2str(nIt)])
  %
  %Plot the spectrum close to the ML estimate
  %
  f = fRef+(-10:.05:10)/st.K;
  sp = nufft(st,f);
  plot(f,20*log10(abs(sp)/M))

  hold on
  stem(f1,20*log10(A1),'r')
  hold off

  grid on
  set(gca,'YLim',[-40,10])
  xlabel('f')
  ylabel('Amplitude (dB)')
  text(fRef+2/st.K,5,{['Freq. ' num2str(f1)],['Ampl. ' num2str(A1)],...
      ['Phase. (rad) ' num2str(Phi1)]})
  %

  R = input('Press Enter for a new trial, q to quit: ','s');

  if any(strcmp(R,{'q','Q'}))
    break
  end

end

%------------------------------------------------------
%------------------------------------------------------

function [f1,A1,Phi1,L1,Deg,nIt] = MLEstimate(st)

%This function finds out the peak of a DFT. Its input is the struct produced by the nufft function. 

%Find the spectral peak in the frequency grid.
[pr11,ind] = max(abs(st.aF).^2);

%Compute its frequency.
f0 = (ind-1)/st.K;
f0 = f0+ceil(-0.5-f0);

%L0 contains the differentials of the ML cost function of orders 0,1, 2 ...
%at f0. 
L0 = LDiff(nufft(st,f0,2));

nItMax = 20;
err = 10.^-12;
nIt = 0;

%Call Iter (Newton iteration) until there is no change or the solution is degenerate.
while 1
  [Deg,pr11,f1,L1] = Iter(st,f0,L0);
  
  if Deg || abs(f1-f0) <= err || nIt >= nItMax
    break
  end
  
  f0 = f1;
  L0 = L1;
  nIt = nIt+1;
end

pr = nufft(st,f1);
A1 = abs(pr)/st.M;
Phi1 = angle(pr); 

function L = LDiff(c)

%Computes the differentials of the ML cost function from the nufft differentials.

L = [real(c(1)*conj(c(1))); 2*real(conj(c(1))*c(2)); ...
  2*real(conj(c(2))*c(2)+conj(c(1))*c(3))];

function [Deg,nCut,f1,L1] = Iter(st,f0,L0)

%Performs the Newton iteration, cutting the Newton direction if necessary.

f1 = f0-L0(2)/L0(3); %Newton iteration.
f1 = f1+ceil(-0.5-f1); %Make it fall in [-0.5,0.5].
L1 = LDiff(nufft(st,f1,2)); %Compute the differentials at f1.

nCutMax = 20; %The difference f1-f0 will be reduced in the sequel if there is no ...
              %improvement in the cost function
nCut = 0;
d = 1;

%In this loop the difference f1-f0 is halved a maximum of nCutMax times until there is ...
%an increase in the cost function

while L1(1) < L0(1) && nCut < nCutMax
  d = d/2;
  f1 = f0-d*L0(2)/L0(3);
  f1 = f1+ceil(-0.5-f1);
  L1 = LDiff(nufft(st,f1,2)); 
  nCut = nCut+1;
end

Deg = nCut >= nCutMax;

%----------------------------------------------------------------
%----------------------------------------------------------------

%Non-uniform FFT function. This code is discussed in 
%   http://www.dsprelated.com/showcode/159.php

function varargout = nufft(varargin)

%Author: J. Selva
%Date: April 2011.
%
%nufft computes the FFT at arbitrary frequencies using barycentric interpolation.
%
%For the interpolation method, see
%
%J. Selva, 'Design of barycentric interpolator for uniform and nonuniform sampling
%  grids', IEEE Trans. on Signal Processing, vol. 58, n. 3, pp. 1618-1627, March 2010.
%
%Usage:
%
% -First call nufft with the vector of samples as argument,
%
%       st = nufft(a); %a is the vector of samples.
%
%    The output is an struct. The field st.aF is the DFT of vector a, sampled with spacing
%    1/st.K. 
%
%    The DFT is defined using
%
%        A_k = \sum_{n=m_1}^{m_1+M-1} a_n e^{-j 2 \pi n f}
%
%    where m_1=-ceil(M/2) and M the number of samples. 
%
% -Then call nufft with sintax 
%
%        out = nufft(st,f,nd,method)
%
%   where
%
%     f: Vector of frequencies where the DFT is evaluated. Its elements must follow
%        abs(f(k))<=0.5
%
%     nd: Derivative order. nufft computes derivatives up to order nd. At the output 
%         out(p,q) is the (q-1)-order derivative of the DFT at frequency f(p). The first
%         column of out contains the zero-order derivatives, i.e, the values of the DFT at
%         frequencies in vector f. 
%     method: If omitted, it is 'baryVec'. Four methods have been implemented:
%
%         +'baryLoop': The output is interpolated using the barycentric method and a loop
%           implementation.
%
%         +'baryVec': The output is interpolated using the barycentric method and a
%          vectorized implementation.
%
%         +'directLoop': The output is computed using the DFT sum directly and a loop
%           implementation.
%  
%         +'directVec': The output is computed using the DFT sum directly and a vectorized
%           implementation.

 
if ~isstruct(varargin{1})
  st = [];
  st.a = varargin{1};
  st.a = st.a(:);
  
  st.M = length(st.a);
  st.m1 = -ceil(st.M/2);
  st.K =  2^ceil(log2(st.M)+1);
  st.aF = fft([st.a(-st.m1+1:st.M); zeros(st.K-st.M,1); st.a(1:-st.m1)]);
  
  errb = 10.^-13; %This is the interpolation accuracy. You can set any larger value, and ...
      %this would reduce st.P. The number of interpolation samples is 2*st.P+1.
  
  st.T = 1/st.K;
  st.B = -2*st.m1;
  st.P = ceil(acsch(errb)/(pi*(1-st.B*st.T)));
  st.vt = MidElementToEnd((-st.P:st.P)*st.T);
  
  st.alpha =  MidElementToEnd(baryWeights(st.T,st.B,st.P));
  
  varargout = {st};
  
  return
end

[st,f] = deal(varargin{1:2});
  
nd = 0;
   
if nargin > 2
  nd = varargin{3};
end
  
method = 'baryVec';
  
if nargin > 3
  method = varargin{4};
end
  
Nf = length(f);
out = zeros(Nf,nd+1);

switch method
    
  case 'baryLoop' %Interpolated method using loops

    for k = 1:length(f)

      x = f(k);
      
      n = floor(x/st.T+0.5);
      u = x - n * st.T;

      da = MidElementToEnd(st.aF(1+mod(n-st.P:n+st.P,st.K)).');

      out(k,:) = DerBarycentricInterp3(st.alpha,da,st.vt,u,nd);

    end

  case 'baryVec' %Vectorized interpolated method
    
    f = f(:);
    Nf = length(f);
    
    n = floor(f/st.T+0.5);
    u = f - n * st.T;
    
    pr = [-st.P:-1 , 1:st.P , 0];
    
    ms = st.aF(1+mod(n(:,ones(1,2*st.P+1)) + pr(ones(Nf,1),:),st.K));
    
    if length(f) == 1
      ms = ms.';
    end
    
    out = DerBarycentricInterp3Vec(st.alpha,ms,st.vt,u,nd);
    
  case 'directLoop' % Direct method using loops
      
    for k = 1:length(f)
      out(k,1) = 0;
	
      for r = st.m1:st.m1+st.M-1
	out(k,1) = out(k,1)+exp(-1i*2*pi*r*f(k))*st.a(r-st.m1+1);
      end
	
      for kd = 1:nd
	out(k,kd+1) = 0;
	
	for r = st.m1:st.m1+st.M-1
	  out(k,kd+1) = out(k,kd+1) + ...
	      ((-1i*2*pi*r).^kd .* exp(-1i*2*pi*r*f(k)))*st.a(r-st.m1+1);
	end
      
      end
    end  
    
  case 'directVec' %Vectorized direct method
	
    for k = 1:length(f)
      out(k,1) = exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k))*st.a;
      
      for kd = 1:nd
	out(k,kd+1) = ...
	    ((-1i*2*pi*(st.m1:st.m1+st.M-1)).^kd .* ...
	    exp(-1i*2*pi*(st.m1:st.m1+st.M-1)*f(k)))*st.a;
      end
      
    end
   
end
 
varargout = {out};

function v = MidElementToEnd(v)

ind = ceil(length(v)/2);
v = [v(1:ind-1),v(ind+1:end),v(ind)];

function v = APPulse(t,B,TSL)

v = real(sinc(B*sqrt(t.^2-TSL^2)))/real(sinc(1i*pi*B*TSL));

function w = baryWeights(T,B,P)

vt = (-P:P)*T;
g = APPulse(vt,1/T-B,T*(P+1));
gam = gamma(vt/T+P+1) .* gamma(-vt/T+P+1) .* g;

N = length(vt);
LD = ones(1,N);

for k = 1:N
  LD(k) = prod(vt(k)-vt(1:k-1))* prod(vt(k)-vt(k+1:N));
end

w = gam ./ LD;
w = w / max(abs(w));

function out = DerBarycentricInterp3(alpha,s,t,tau,nd)

vD = 0;
Nt = length(t);
LF = [1,zeros(1,Nt-1)];
out = zeros(1,nd+1);

for k = 1:Nt-1
  vD = vD*(tau-t(k))+alpha(k)*LF(k);
  LF(k+1) = LF(k)*(tau-t(k));
end

vD = vD * (tau-t(Nt)) + alpha(Nt) * LF(Nt);

z = s;

for kd = 0:nd

  z1 = z-z(end); cN = 0;

  for k = 1:Nt-1
    cN = cN * (tau-t(k))+z1(k)*alpha(k)*LF(k);
  end
  cN = cN/vD;
  
  ztau = z(end)+(tau-t(end))*cN;
  out(kd+1) = ztau;
  
  if kd < nd
    z = [ (z(1:end-1)-ztau)./(t(1:end-1)-tau) , cN ] * (kd+1);
  end
  
end

function out = DerBarycentricInterp3Vec(alpha,zL,t,tauL,nd)

NtauL = size(tauL,1);
LBtau = 200; %Block size for vectorized processing. Change this to adjust the memory 
             %usage. 

NBlock = 1+floor(NtauL/LBtau);

Nt = length(t);

out = zeros(NtauL,nd+1);

for r = 0:NBlock-1
  ind1 = 1+r*LBtau;
  ind2 = min([(r+1)*LBtau,NtauL]);
  
  Ntau = ind2-ind1+1;
  z = zL(ind1:ind2,:);
  tau = tauL(ind1:ind2);
  
  vD = zeros(Ntau,1);
  
  LF = [1,zeros(1,Nt-1)];
  LF = LF(ones(Ntau,1),:);

  for k = 1:Nt-1
    vD = vD .* (tau-t(k))+alpha(k)*LF(:,k);
    LF(:,k+1) = LF(:,k) .* (tau-t(k));
  end

  vD = vD .* (tau-t(Nt)) + alpha(Nt) * LF(:,Nt);

  for kd = 0:nd

    pr = z(:,end); z1 = z-pr(:,ones(1,Nt)); cN = zeros(Ntau,1);

    for k = 1:Nt-1
      cN = cN .* (tau-t(k)) + z1(:,k) .* alpha(k) .* LF(:,k);
    end
    cN = cN ./ vD;
  
    ztau = z(:,end)+(tau-t(end)) .* cN;
    out(ind1:ind2,kd+1) = ztau;
  
    if kd < nd
      pr1 = ztau;
      pr1 = z(:,1:end-1) - pr1(:,ones(1,Nt-1));
      
      pr2 = t(1:end-1);
      pr2 = pr2(ones(Ntau,1),:)-tau(:,ones(1,Nt-1));
      
      z = [pr1 ./ pr2, cN] * (kd+1);
      
    end
  end
  
end

Discrete Cosine Transform 1D Matrix Processing

Emmanuel ● April 16, 2011 ● Coded in Matlab 
function [Sdct] = matrix_dct(P)    %where P is the vector to apply the dct
N=length(P);
Cf=[1/sqrt(2),ones(1,N-1)];
S=zeros(1,N);
for f=0: N-1;
    for x=0: N-1;
        W(f+1, x+1)=cos((2*x+1)*pi*f/(2*N));  %matrix kernel
    end
end
Sdct=W*P';   
Sdct=sqrt(2/N)*Sdct.*Cf';     %constant product

Discrete Cosine Transform 2D Serial Processing

Emmanuel ● April 16, 2011 ● Coded in Matlab 
function [S] = dct_2d(P)    %where P is the matrix to apply the dct
N=length(P);
Cu=[1/sqrt(2),ones(1,N-1)];   %Constant coeficients d1
Cv=[1/sqrt(2),ones(1,N-1)];   %Constant coeficients d2
S=zeros(N,N);
for v=0: N-1;
    for y=0: N-1;
        for u=0: N-1;
            for x=0: N-1;    %serial processing
                S(u+1,v+1)=S(u+1,v+1)+Cu(u+1)*Cv(v+1)*P(x+1,y+1)*cos((2*x+1)*pi*u/(2*N))*cos((2*y+1)*pi*v/(2*N));
            end;
        end;
    end;
end;
S=(2/N)*S;  %final result

Frequency Response Creator

Ron ● April 11, 2011 ● Coded in Matlab 
[header, s] = hdrload('data_file.txt');

num = 2;
x = zeros(ceil(44100/num),1);                 %channel 1
y = zeros(ceil(44100/num),1);                 %channel 2
 
Fs = 44100;                     % Sampling frequency
T = 1/Fs;                       % Sample time
L = Fs/num;                  	  % Length of window
g = floor(num*length(s)/Fs);    % number of windows
H = 0;
count = 50;
max_index = 1;                  % pointer for max value
max = 0;                        % max value
max_index2 = 1;                 % pointer for max value
max2 = 0;                       % max value
 
z = zeros(g,2);
z2 = zeros(g,2);
 
for i = 1:1:g-3
    
% creating windows for channel one and two
    for j = 1:1:L
       k = j + L*i;
       x(j,1) = s(k+i,1);
       y(j,1) = s(k+i,2);
    end
 
NFFT = 2^nextpow2(L);            % Next power of 2 from length of y
YY = fft(y,NFFT)/L;
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1); % creating frequency vector
Y2 = 2*abs(Y(1:NFFT/2+1));
YY2 = 2*abs(YY(1:NFFT/2+1));
 
%%channel 1
max = Y(max_index,1);
 
%%channel 2 
if(max_index2 < length(YY2))
max2 = YY(max_index2,1);
 
b = length(YY2) - max_index2 - 2*count;
if b > 0
    a = max_index2 + 2*count; 
else
    a = length(YY2);
end
 
for j = max_index2:1:(a-1)
    if max2 < YY2(j+1,1)
        max2 = YY2(j+1,1);
        max = Y2(j+1,1);
        max_index2 = j+1;
    end
end
end
    z2(i+1,1) = abs(max2);
    z2(i+1,2) = f(1,max_index2);
    z(i+1,1) = abs(max);
    z(i+1,2) = f(1,max_index2);
    if(max_index2 > 4*count)
    max_index2 = max_index2 - count;
    end

% RECORD MOVIE
% semilogy(z2(:,2),z2(:,1),f,YY2,z(:,2),z(:,1),f,Y2);
% axis([0 20000 10^-3 Inf]);

end
 
semilogy(z2(:,2),z2(:,1),z(:,2),z(:,1));
axis([0 20000 0.00001 Inf]);

OFDM symbol generator

kaz - ● April 4, 2011●3 comments ● Coded in Matlab 
%bandwidth = nCar/nFFT *Fs

clear all;
nData = 200000;                %vector length
ModulationType = '64QAM';
nFFT = 2048;                   %fft size
nCar = 1200;                   %number of active carriers
nPrefix = 144;
Fs = 10;                       %sampling rate
       
%compute number of ofdm symbols  
nSymbols = ceil(nData/(nFFT+nPrefix));
   
%generate random data symbols
switch ModulationType                              
    case 'QPSK',  alphabet = [-1,1];
    case '16QAM', alphabet = [-3,-1,1,3];
    case '64QAM', alphabet = [-7,-5,-3,-1,1,3,5,7];
end
   
tx = [];
for k = 1:nSymbols    
   QAM = complex(randsrc(nCar,1,alphabet),randsrc(nCar,1,alphabet));
 
   fft_in = zeros(1,nFFT);
   fft_in(1:nCar/2) = QAM(1:nCar/2);
   fft_in(end-nCar/2+1:end) = QAM(nCar/2+1:end);
   fft_out = ifft(fft_in);

   
   %concatenate ofdm symbol plus its prefix into tx vector
   prefix = fft_out(nFFT-nPrefix+1:nFFT);          
   tx = [tx prefix fft_out];
end

[P, F] = pwelch(tx, hann(2^16), 0, 2^16, Fs);
P = fftshift(P);
F = F - max(F)/2;
plot(F,10*log10(P));