Peak/Notch Filter Design

May 21, 20113 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

May 18, 20111 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

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 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

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;

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()) {
/*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

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 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

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[.

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)],...
%

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

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

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

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

April 4, 20113 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));