% *******************************************************
% delay-matching between two signals (complex/real-valued)
% M. Nentwig
%
% * matches the continuous-time equivalent waveforms
% of the signal vectors (reconstruction at Nyquist limit =>
% ideal lowpass filter)
% * Signals are considered cyclic. Use arbitrary-length
% zero-padding to turn a one-shot signal into a cyclic one.
%
% * output:
% => coeff: complex scaling factor that scales 'ref' into 'signal'
% => delay 'deltaN' in units of samples (subsample resolution)
% apply both to minimize the least-square residual
% => 'shiftedRef': a shifted and scaled version of 'ref' that
% matches 'signal'
% => (signal - shiftedRef) gives the residual (vector error)
%
% *******************************************************
function [coeff, shiftedRef, deltaN] = fitSignal_FFT(signal, ref)
n=length(signal);
% xyz_FD: Frequency Domain
% xyz_TD: Time Domain
% all references to 'time' and 'frequency' are for illustration only
forceReal = isreal(signal) && isreal(ref);
% *******************************************************
% Calculate the frequency that corresponds to each FFT bin
% [-0.5..0.5[
% *******************************************************
binFreq=(mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
% *******************************************************
% Delay calculation starts:
% Convert to frequency domain...
% *******************************************************
sig_FD = fft(signal);
ref_FD = fft(ref, n);
% *******************************************************
% ... calculate crosscorrelation between
% signal and reference...
% *******************************************************
u=sig_FD .* conj(ref_FD);
if mod(n, 2) == 0
% for an even sized FFT the center bin represents a signal
% [-1 1 -1 1 ...] (subject to interpretation). It cannot be delayed.
% The frequency component is therefore excluded from the calculation.
u(length(u)/2+1)=0;
end
Xcor=abs(ifft(u));
% figure(); plot(abs(Xcor));
% *******************************************************
% Each bin in Xcor corresponds to a given delay in samples.
% The bin with the highest absolute value corresponds to
% the delay where maximum correlation occurs.
% *******************************************************
integerDelay = find(Xcor==max(Xcor));
% (1): in case there are several bitwise identical peaks, use the first one
% Minus one: Delay 0 appears in bin 1
integerDelay=integerDelay(1)-1;
% Fourier transform of a pulse shifted by one sample
rotN = exp(2i*pi*integerDelay .* binFreq);
uDelayPhase = -2*pi*binFreq;
% *******************************************************
% Since the signal was multiplied with the conjugate of the
% reference, the phase is rotated back to 0 degrees in case
% of no delay. Delay appears as linear increase in phase, but
% it has discontinuities.
% Use the known phase (with +/- 1/2 sample accuracy) to
% rotate back the phase. This removes the discontinuities.
% *******************************************************
% figure(); plot(angle(u)); title('phase before rotation');
u=u .* rotN;
% figure(); plot(angle(u)); title('phase after rotation');
% *******************************************************
% Obtain the delay using linear least mean squares fit
% The phase is weighted according to the amplitude.
% This suppresses the error caused by frequencies with
% little power, that may have radically different phase.
% *******************************************************
weight = abs(u);
constRotPhase = 1 .* weight;
uDelayPhase = uDelayPhase .* weight;
ang = angle(u) .* weight;
r = [constRotPhase; uDelayPhase] .' \ ang.'; %linear mean square
%rotPhase=r(1); % constant phase rotation, not used.
% the same will be obtained via the phase of 'coeff' further down
fractionalDelay=r(2);
% *******************************************************
% Finally, the total delay is the sum of integer part and
% fractional part.
% *******************************************************
deltaN = integerDelay + fractionalDelay;
% *******************************************************
% provide shifted and scaled 'ref' signal
% *******************************************************
% this is effectively time-convolution with a unit pulse shifted by deltaN
rotN = exp(-2i*pi*deltaN .* binFreq);
ref_FD = ref_FD .* rotN;
shiftedRef = ifft(ref_FD);
% *******************************************************
% Again, crosscorrelation with the now time-aligned signal
% *******************************************************
coeff=sum(signal .* conj(shiftedRef)) / sum(shiftedRef .* conj(shiftedRef));
shiftedRef=shiftedRef * coeff;
if forceReal
shiftedRef = real(shiftedRef);
end
end
#
# rfft.py
#
# This file contains a recursive version of the fast-fourier transform and
# support test functions. This module utilizes the numpy (numpy.scipy.org)
# library.
#
# References
# - http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
# - "A Simple and Efficient FFT Implementation in C++", by Vlodymyr Myrnyy
import numpy
from numpy.fft import fft
from numpy import sin, cos, pi, ones, zeros, arange, r_, sqrt, mean
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def rFFT(x):
"""
Recursive FFT implementation.
References
-- http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
-- "A Simple and Efficient FFT Implementation in C++"
by Vlodymyr Myrnyy
"""
n = len(x)
if (n == 1):
return x
w = getTwiddle(n)
m = n/2;
X = ones(m, float)*1j
Y = ones(m, float)*1j
for k in range(m):
X[k] = x[2*k]
Y[k] = x[2*k + 1]
X = rFFT(X)
Y = rFFT(Y)
F = ones(n, float)*1j
for k in range(n):
i = (k%m)
F[k] = X[i] + w[k] * Y[i]
return F
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def getTwiddle(NFFT=8):
"""Generate the twiddle factors"""
W = r_[[1.0 + 1.0j]*NFFT]
for k in range(NFFT):
W[k] = cos(2.0*pi*k/NFFT) - 1.0j*sin(2.0*pi*k/NFFT)
return W
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def DFT(x, N=8):
"""
Use the direct definition of DFT for verification
"""
y = [1.0 + 1.0j]*N
y = r_[y]
for n in range(N):
wsum = 0 + 0j;
for k in range(N):
wsum = wsum + (cos(2*pi*k*n/N) - (1.0j * sin(2*pi*k*n/N)))*x[k]
y[n] = wsum
return y
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def test_rfft(N = 64, # FFT order to test
nStart = 0.2, # Note aliased signal included
nStep = 2.1, # Samples per period step
pStep = pi/4, # Phase step size
limErr = 10e-12, # Error limit to check
maxErr = 0 # Max difference
):
"""
Use the built in numpy FFT functions and the direct
implemenation of the DFT to verify the recursive FFT.
This testbench verifies the different implementations are within
a certain limit. Because of the different implemenations the values
could be slightly off (computer representation calculation error).
"""
# Use test signal nStart:nStep:N samples per cycle
for s in arange(nStart, N+nStep, nStep):
for p in arange(0, pi+pStep, pStep):
n = arange(N, 0, -1)
x = cos(2*pi*n/s + p)
xDFT = DFT(x,N)
nFFT = fft(x,N)
xFFT = rFFT(x)
rmsErrD = sqrt(mean(abs(xDFT - xFFT))**2)
rmsErrN = sqrt(mean(abs(nFFT - xFFT))**2)
if rmsErrD > limErr or rmsErrN > limErr:
print s, p, "Error!", rmsErrD, rmsErrN
print xDFT
print nFFT
print xFFT
if rmsErrD > maxErr:
maxErr = rmsErrD
elif rmsErrN > maxErr:
maxErr = rmsErrN
print "N %d maxErr = %f " % (N,maxErr)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If the module is run test a bunch of different size FFTs
if __name__ == '__main__':
# The following is fairly exhaustive and will take some time
# to run.
tv = 2**arange(1,12)
for nfft in tv:
test_rfft(N=nfft)
% 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)';
/* -------------- MovAvg.h begins -------------- */
#ifndef __MOVAVG_H__
#define __MOVAVG_H__
#include "types.h"
typedef struct
{
unsigned int decimation;
unsigned int deccnt;
unsigned int BufHeadPtr;
unsigned int BufLength;
Frac16 multiplier;
Frac16 ykm1;
}MovAvgparams;
/*---------------------------------------------------------------------------;
; Initialize buffer length, pointer and decimation divider of moving average ;
; filter; initialize memory buffer of moving average filter to all zeroes. ;
; ;
; Input: Y0 = buffer length ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = buffer length (filter order) ;
; R3 = pointer to storage buffer for the past input samples ;
; ;
; Output: None ;
; ;
; Registers modified: R2 ;
;---------------------------------------------------------------------------*/
asm void MovAvgInit(unsigned int length, MovAvgparams *params, Frac16 *data);
/*---------------------------------------------------------------------------;
; This function implements a moving average Nth order FIR filter. ;
; The result yk is computed as the sum of the current and the past N input ;
; samples multiplied by the multiplier value. Canonical moving average ;
; filtering is obtained if multiplier = 1/(N+1). ;
; For optimal results it is advised to call this function with saturation ;
; disabled. ;
; ;
; Input: Y0 = xk (input sample) ;
; (R2) = decimation (1 = no decimation, 2 = discard every other ;
; sample, and so on) ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = N (length of buffer and filter order) ;
; (R2+4) = fractional multiplier ;
; (R2+5) = yk-1 (past output) ;
; R3 = Pointer to storage buffer for the N past input samples ;
; ;
; Output: Y0 = yk ;
; ;
; Registers modified: A, X0, Y, R0, R3 ;
;---------------------------------------------------------------------------*/
asm Frac16 MovAvg(Frac16 xk, MovAvgparams *params, Frac16 *data);
#endif //ifndef __MOVAVG_H__
/* -------------- MovAvg.h ends -------------- */
/* -------------- MovAvg.c begins -------------- */
/*---------------------------------------------------------------------------;
; Initialize buffer length, pointer and decimation divider of moving average ;
; filter; initialize memory buffer of moving average filter to all zeroes. ;
; ;
; Input: Y0 = buffer length ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = buffer length (filter order) ;
; R3 = pointer to storage buffer for the past input samples ;
; ;
; Output: None ;
; ;
; Registers modified: R2 ;
;---------------------------------------------------------------------------*/
asm void MovAvgInit(unsigned int length, MovAvgparams *params, Frac16 *data){
MOVE.W #0001,X:(R2+1) //Initialize decimation divider
CLR.W X:(R2+2) //Initialize buffer head pointer
MOVE.W Y0,X:(R2+3) //Initialize buffer length
REP Y0 //Initialize past data to zero
CLR.W X:(R3)+
RTS //Return from subroutine
}
/*---------------------------------------------------------------------------;
; This function implements a moving average Nth order FIR filter. ;
; The result yk is computed as the sum of the current and the past N input ;
; samples multiplied by the multiplier value. Canonical moving average ;
; filtering is obtained if multiplier = 1/(N+1). ;
; For optimal results it is advised to call this function with saturation ;
; disabled. ;
; ;
; Input: Y0 = xk (input sample) ;
; (R2) = decimation (1 = no decimation, 2 = discard every other ;
; sample, and so on) ;
; (R2+1) = decimation counter (used internally for data decimation) ;
; (R2+2) = index of buffer head ;
; (R2+3) = N (length of buffer and filter order) ;
; (R2+4) = fractional multiplier ;
; (R2+5) = yk-1 (past output) ;
; R3 = Pointer to storage buffer for the N past input samples ;
; ;
; Output: Y0 = yk ;
; ;
; Registers modified: A, X0, Y, R2, R3 ;
;---------------------------------------------------------------------------*/
asm Frac16 MovAvg(Frac16 xk, MovAvgparams *params, Frac16 *data){
DEC.W X:(R2+1) //Decrement decimation counter
BEQ DoMovAvg //Decremented to zero?
MOVE.W X:(R2+5),Y0 //No, output old value and exit
BRA MovAvgDone
DoMovAvg: MOVE.W X:(R2)+,X0 //Yes, reload decimation counter
MOVE.W X0,X:(R2)+
ADDA #1,SP //Preserve N (needed if called
MOVE.W N,X:(SP) //from an ISR)
MOVE.W X:(R2)+,A //Load index of buffer head into
MOVE.W A1,N //address modifier register
INC.W A X:(R2)+,Y1 //Point to next buffer location
//and load buffer length
CMP.W Y1,A //If buffer head index = buffer
BNE MAIndexOK //length, then reset buffer head
CLR.W A1 //index
MAIndexOK: MOVE.W A1,X:(R2-2) //Save buffer head index
MOVE.W X:(R3+N),A1 //Load oldest buffer sample
MOVE.W Y0,X:(R3+N) //Save new input sample
MOVE.W X:(R2)+,Y0 //Load multiplier
MPY A1,Y0,A X:(R3)+,X0 //Accumulate oldest sample
REP Y1 //Accumulate the latest N samples
MAC X0,Y0,A X:(R3)+,X0
MOVE.W A,Y0 //Return accumulation result
MOVE.W Y0,X:(R2) //Save new output sample
MOVE.W X:(SP)-,X0 //Restore N
MOVE.W X0,N
MovAvgDone: RTS //Return from subroutine
}
/* -------------- MovAvg.c ends -------------- */
/* -------------- Usage example begins ------------- */
/* application specific constants */
#define MovAvgOrder 40 //Order of moving average filter
#define MovAvgDecimation 100 //Decimation factor for m.a.
/* application specific includes */
#include "MovAvg.h"
/* global variables */
MovAvgparams MAparams;
Frac16 MAdata[MovAvgOrder];
int FiltIn,MAOut;
/* initializations */
MAparams.decimation=MovAvgDecimation;
MAparams.multiplier=FRAC16((float)1/(MovAvgOrder+1));
MovAvgInit(MovAvgOrder,&MAparams,(Frac16*)&MAdata); //Initialize m.a. filter
/* Moving average filter function call */
MAOut=MovAvg(FiltIn,&MAparams,(Frac16*)&MAdata); //Call moving av. filter
/* -------------- Usage example ends ------------- */
/* David Valencia at DSPRelated.com
This example shows how to use the External Memory Interface
expansion port of Spectrum Digital DSK6713 card as
a set of digital imputs
The use of buffer ICs is mandatory, as these pins can't deliver
too much current */
#define OUTPUT 0xA0000000 // EMIF ADDRESS (updated 3/11/2010)
int *output = (int*)OUTPUT; // Declare a pointer to EMIF's address
void main()
{
// To see how the value is written to the EMIF pins,
// you may prefer to run the program step-by-step
// Write a 32-bit digital value to the pins
*output = 0x48120A00; // Test value
}
function [D1,D2] = Dec_OptimTwoStage_Factors(D_Total,Passband_width,Fstop)
%
% When breaking a single large decimation factor (D_Total) into
% two stages of decimation (to minimize the orders of the
% necessary lowpass digital filtering), an optimum first
% stage of decimation factor (D1) can be found. The second
% stage's decimation factor is then D_Total/D1.
%
% Inputs:
% D_Total = original single-stage decimation factor.
% Passband_width = desired passband width of original
% single-stage lowpass filter (Hz).
% Fstop = desired beginning of stopband freq of original
% single-stage lowpass filter (Hz).
%
% Outputs:
% D1 = optimum first-stage decimation factor.
% D2 = optimum second-stage decimation factor.
%
% Example: Assume you want to decimate a sequence by D_Total=100,
% your original single-stage lowpass filter's passband
% width and stopband frequency are 1800 and 2200 Hz
% respectively. We use:
%
% [D1,D2] = Dec_OptimTwoStage_Factors(100, 1800, 2200)
%
% giving us the desired D1 = 25, and D2 = 4.
% (That is, D1*D2 = 25*4 = 100 = D_Total.)
%
% Author: Richard Lyons, November, 2011
% Start of processing
DeltaF = (Fstop -Passband_width)/Fstop;
Radical = (D_Total*DeltaF)/(2 -DeltaF); % Terms under sqrt sign.
Numer = 2*D_Total*(1 -sqrt(Radical)); % Numerator.
Denom = 2 -DeltaF*(D_Total + 1); % Denominator.
D1_estimated = Numer/Denom; %Optimum D1 factor, but not
% an integer.
% Find all factors of 'D_Total'
Factors_of_D_Total = Find_Factors_of_D_Total(D_Total);
% Find the one factor of 'D_Total' that's closest to 'D1_estimated'
Temp = abs(Factors_of_D_Total -D1_estimated);
Index_of_Min = find(Temp == min(Temp)); % Index of optim. D1
D1 = Factors_of_D_Total(Index_of_Min); % Optimum first
% decimation factor
D2 = D_Total/D1; % Second decimation factor.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allfactors] = Find_Factors_of_D_Total(X)
% Finds all integer factors of intger X.
% Filename was originally 'listfactors.m', written
% by Jeff Miller. Downloaded from Matlab Central.
[b,m,n] = unique(factor(X));
%b is all prime factors
occurences = [m(1) diff(m)];
current_factors = [b(1).^(0:occurences(1))]';
for index = 2:length(b)
currentprime = b(index).^(0:occurences(index));
current_factors = current_factors*currentprime;
current_factors = current_factors(:);
end
allfactors = sort(current_factors);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//Hamming Encoding
//H(7,4)
//Code Word Length = 7, Message Word length = 4, Parity bits =3
//clear;
close;
clc;
//Getting Message Word
m3 = input('Enter the 1 bit(MSb) of message word');
m2 = input('Enter the 2 bit of message word');
m1 = input('Enter the 3 bit of message word');
m0 = input('Enter the 4 bit(LSb) of message word');
//Generating Parity bits
for i = 1:(2^4)
b2(i) = xor(m0(i),xor(m3(i),m1(i)));
b1(i) = xor(m1(i),xor(m2(i),m3(i)));
b0(i) = xor(m0(i),xor(m1(i),m2(i)));
m(i,:) = [m3(i) m2(i) m1(i) m0(i)];
b(i,:) = [b2(i) b1(i) b0(i)];
end
C = [b m];
disp('___________________________________________________________')
for i = 1:2^4
disp(i)
disp(m(i,:),'Message Word')
disp(b(i,:),'Parity Bits')
disp(C(i,:),'CodeWord')
disp(" ");
disp(" ");
end
disp('___________________________________________________________')
//disp(m b C)
//Result
//Enter the 1 bit(MSb) of message word [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1];
//Enter the 2 bit of message word [0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1];
//Enter the 3 bit of message word [0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1];
//Enter the 4 bit(LSb) of message word [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1];
//Determination of spectrum of a signal
//With maximum normalized frequency f = 0.4
//using Blackmann window
clear all;
close;
clc;
N = 11;
cfreq = [0.4 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
wft; // Time domain filter coefficients
wfm; // Frequency domain filter values
fr; // Frequency sample points
for n = 1:N
h_balckmann(n)=0.42-0.5*cos(2*%pi*n/(N-1))+0.08*cos(4*%pi*n/(N-1));
wft_blmn(n) = wft(n)*h_balckmann(n);
end
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
plot2d(fr,WFM_blmn_dB)
xtitle('Frequency Response of Blackmann window Filtered output N = 11','Frequency in cycles per samples f','Energy density in dB')
//Evaluating power spectrum of a discrete sequence
//Using N-point DFT
clear all;
clc;
close;
N =32; //Number of samples in given sequence
n =0:N-1;
delta_f = [0.06,0.01];//frequency separation
x1 = sin(2*%pi*0.315*n)+cos(2*%pi*(0.315+delta_f(1))*n);
x2 = sin(2*%pi*0.315*n)+cos(2*%pi*(0.315+delta_f(2))*n);
L = [8,64,256,512];
k1 = 0:L(1)-1;
k2 = 0:L(2)-1;
k3 = 0:L(3)-1;
k4 = 0:L(4)-1;
fk1 = k1./L(1);
fk2 = k2./L(2);
fk3 = k3./L(3);
fk4 = k4./L(4);
for i =1:length(fk1)
Pxx1_fk1(i) = 0;
Pxx2_fk1(i) = 0;
for m = 1:N
Pxx1_fk1(i)=Pxx1_fk1(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk1(i));
Pxx2_fk1(i)=Pxx1_fk1(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk1(i));
end
Pxx1_fk1(i) = (Pxx1_fk1(i)^2)/N;
Pxx2_fk1(i) = (Pxx2_fk1(i)^2)/N;
end
for i =1:length(fk2)
Pxx1_fk2(i) = 0;
Pxx2_fk2(i) = 0;
for m = 1:N
Pxx1_fk2(i)=Pxx1_fk2(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk2(i));
Pxx2_fk2(i)=Pxx1_fk2(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk2(i));
end
Pxx1_fk2(i) = (Pxx1_fk2(i)^2)/N;
Pxx2_fk2(i) = (Pxx1_fk2(i)^2)/N;
end
for i =1:length(fk3)
Pxx1_fk3(i) = 0;
Pxx2_fk3(i) = 0;
for m = 1:N
Pxx1_fk3(i) =Pxx1_fk3(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk3(i));
Pxx2_fk3(i) =Pxx1_fk3(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk3(i));
end
Pxx1_fk3(i) = (Pxx1_fk3(i)^2)/N;
Pxx2_fk3(i) = (Pxx1_fk3(i)^2)/N;
end
for i =1:length(fk4)
Pxx1_fk4(i) = 0;
Pxx2_fk4(i) = 0;
for m = 1:N
Pxx1_fk4(i) =Pxx1_fk4(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk4(i));
Pxx2_fk4(i) =Pxx1_fk4(i)+x1(m)*exp(-sqrt(-1)*2*%pi*(m-1)*fk4(i));
end
Pxx1_fk4(i) = (Pxx1_fk4(i)^2)/N;
Pxx2_fk4(i) = (Pxx1_fk4(i)^2)/N;
end
figure
title('for frequency separation = 0.06')
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx1_fk1))
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx1_fk2))
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx1_fk3))
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx1_fk4))
figure
title('for frequency separation = 0.01')
subplot(2,2,1)
plot2d3('gnn',k1,abs(Pxx2_fk1))
subplot(2,2,2)
plot2d3('gnn',k2,abs(Pxx2_fk2))
subplot(2,2,3)
plot2d3('gnn',k3,abs(Pxx2_fk3))
subplot(2,2,4)
plot2d3('gnn',k4,abs(Pxx2_fk4))
//Using Digital Filter Transformation, the First order
//Analog IIR Butterworth LPF converted into Digital
//Butterworth HPF
clear all;
clc;
close;
s = poly(0,'s');
Omegac = 0.2*%pi;
H = Omegac/(s+Omegac);
T =1;//Sampling period T = 1 Second
z = poly(0,'z');
Hz_LPF = horner(H,(2/T)*((z-1)/(z+1)));
alpha = -(cos((Omegac+Omegac)/2))/(cos((Omegac-Omegac)/2));
HZ_HPF=horner(Hz_LPF,-(z+alpha)/(1+alpha*z))
HW =frmag(HZ_HPF(2),HZ_HPF(3),512);
W = 0:%pi/511:%pi;
plot(W/%pi,HW)
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(1)
xtitle('Magnitude Response of Single pole HPF Filter Cutoff frequency = 0.2*pi','Normalized Digital Frequency W/pi--->','Magnitude');
%Normalized Least Mean Square Adaptive Filter
%for Echo Cancellation
function [e,h,t] = adapt_filt_nlms(x,B,h,delta,l,l1)
%x = the signal from the speaker
%B = signal through echo path
%h = impulse response of adaptive filter
%l = length of the signal x
%l1 = length of the adaptive filter
%delta = step size
%e = error
%h = adaptive filter impulse response
tic;
for k = 1:150
for n= 1:l
xcap = x((n+l1-1):-1:(n+l1-1)-l1+1);
yout(n) = h*xcap';
e(n) = B(n) - yout(n);
xnorm = 0.000001+(xcap*xcap');
h = h+((delta*e(n))*(xcap/xnorm));
end
eold = 0.0;
for i = 1:l
mesquare = eold+(e(i)^2);
eold = mesquare;
end
if mesquare <=0.0001
break;
end
end
t = toc; %to get the time elapsed
