/* 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
}
# The following is a Python/scipy snippet to generate the
# coefficients for a halfband filter. A halfband filter
# is a filter where the cutoff frequency is Fs/4 and every
# other coeffecient is zero except the cetner tap.
# Note: every other (even except 0) is 0, most of the coefficients
# will be close to zero, force to zero actual
import numpy
from numpy import log10, abs, pi
import scipy
from scipy import signal
import matplotlib
import matplotlib.pyplot
import matplotlib as mpl
# ~~[Filter Design with Parks-McClellan Remez]~~
N = 32 # Filter order
# Filter symetric around 0.25 (where .5 is pi or Fs/2)
bands = numpy.array([0., .22, .28, .5])
h = signal.remez(N+1, bands, [1,0], [1,1])
h[abs(h) <= 1e-4] = 0.
(w,H) = signal.freqz(h)
# ~~[Filter Design with Windowed freq]~~
b = signal.firwin(N+1, 0.5)
b[abs(h) <= 1e-4] = 0.
(wb, Hb) = signal.freqz(b)
# Dump the coefficients for comparison and verification
print(' remez firwin')
print('------------------------------------')
for ii in range(N+1):
print(' tap %2d %-3.6f %-3.6f' % (ii, h[ii], b[ii]))
## ~~[Plotting]~~
# Note: the pylab functions can be used to create plots,
# and these might be easier for beginners or more familiar
# for Matlab users. pylab is a wrapper around lower-level
# MPL artist (pyplot) functions.
fig = mpl.pyplot.figure()
ax0 = fig.add_subplot(211)
ax0.stem(numpy.arange(len(h)), h)
ax0.grid(True)
ax0.set_title('Parks-McClellan (remez) Impulse Response')
ax1 = fig.add_subplot(212)
ax1.stem(numpy.arange(len(b)), b)
ax1.set_title('Windowed Frequency Sampling (firwin) Impulse Response')
ax1.grid(True)
fig.savefig('hb_imp.png')
fig = mpl.pyplot.figure()
ax1 = fig.add_subplot(111)
ax1.plot(w, 20*log10(abs(H)))
ax1.plot(w, 20*log10(abs(Hb)))
ax1.legend(['remez', 'firwin'])
bx = bands*2*pi
ax1.axvspan(bx[1], bx[2], facecolor='0.5', alpha='0.33')
ax1.plot(pi/2, -6, 'go')
ax1.axvline(pi/2, color='g', linestyle='--')
ax1.axis([0,pi,-64,3])
ax1.grid('on')
ax1.set_ylabel('Magnitude (dB)')
ax1.set_xlabel('Normalized Frequency (radians)')
ax1.set_title('Half Band Filter Frequency Response')
fig.savefig('hb_rsp.png')
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ****************************************************************
% RMS phase error from phase noise spectrum
% Markus Nentwig, 2011
%
% - integrates RMS phase error from a phase noise power spectrum
% - generates an example signal and determines the RMS-average
% phase error (alternative method)
% ****************************************************************
function pn_snippet()
close all;
% ****************************************************************
% Phase noise spectrum model
% --------------------------
% Notes:
% - Generally, a phase noise spectrum tends to follow
% |H(f)| = c0 + c1 f^n1 + c2 f^n2 + c3 f^n3 + ...
% Therefore, linear interpolation should be performed in dB
% on a LOGARITHMIC frequency axis.
% - Single-/double sideband definition:
% The phase noise model is defined for -inf <= f <= inf
% in other words, it contributes the given noise density at
% both positive AND negative frequencies.
% Assuming a symmetrical PN spectrum, the model is evaluated
% for |f| => no need to explicitly write out the negative frequency
% side.
% ****************************************************************
% PN model
% first column:
% frequency offset
% second column:
% spectral phase noise density, dB relative to carrier in a 1 Hz
% observation bandwidth
d = [0 -80
10e3 -80
1e6 -140
9e9 -140];
f_Hz = d(:, 1) .';
g_dBc1Hz = d(:, 2) .' -3;
% get RMS phase error by integrating the power spectrum
% (alternative 1)
e1_degRMS = integratePN_degRMS(f_Hz, g_dBc1Hz)
% get RMS phase error based on a test signal
% (alternative 2)
n = 2 ^ 20;
deltaF_Hz = 2;
e2_degRMS = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz)
end
% ****************************************************************
% Integrate the RMS phase error from the power spectrum
% ****************************************************************
function r = integratePN_degRMS(f1_Hz, g1_dBc1Hz)
1;
% integration step size
deltaF_Hz = 100;
% list of integration frequencies
f_Hz = deltaF_Hz:deltaF_Hz:5e6;
% interpolate spectrum on logarithmic frequency, dB scale
% unit is dBc in 1 Hz, relative to a unity carrier
fr_dB = interp1(log(f1_Hz+eps), g1_dBc1Hz, log(f_Hz+eps), 'linear');
% convert to power in 1 Hz, relative to a unity carrier
fr_pwr = 10 .^ (fr_dB/10);
% scale to integration step size
% unit: power in deltaF_Hz, relative to unity carrier
fr_pwr = fr_pwr * deltaF_Hz;
% evaluation frequencies are positive only
% phase noise is two-sided
% (note the IEEE definition: one-half the double sideband PSD)
fr_pwr = fr_pwr * 2;
% sum up relative power over all frequencies
pow_relToCarrier = sum(fr_pwr);
% convert the phase noise power to an RMS magnitude, relative to the carrier
pnMagnitude = sqrt(pow_relToCarrier);
% convert from radians to degrees
r = pnMagnitude * 180 / pi;
end
% ****************************************************************
% Generate a PN signal with the given power spectrum and determine
% the RMS phase error
% ****************************************************************
function r = simulatePN_degRMS(f_Hz, g_dBc1Hz, n, deltaF_Hz);
A = 1; % unity amplitude, arbitrary
% indices of positive-frequency FFT bins.
% Does not include DC
ixPos = 2:(n/2);
% indices of negative-frequency FFT bins.
% Does not include the bin on the Nyquist limit
assert(mod(n, 2) == 0, 'n must be even');
ixNeg = n - ixPos + 2;
% Generate signal in the frequency domain
sig = zeros(1, n);
% positive frequencies
sig(ixPos) = randn(size(ixPos)) + 1i * randn(size(ixPos));
% for purely real-valued signal: conj()
% for purely imaginary-valued signal such as phase noise: -conj()
% Note:
% Small-signals are assumed. For higher "modulation indices",
% phase noise would become a circle in the complex plane
sig(ixNeg) = -conj(sig(ixPos));
% normalize energy to unity amplitude A per bin
sig = sig * A / RMS(sig);
% normalize energy to 0 dBc in 1 Hz
sig = sig * sqrt(deltaF_Hz);
% frequency vector corresponding to the FFT bins (0, positive, negative)
fb_Hz = FFT_freqbase(n, deltaF_Hz);
% interpolate phase noise spectrum on frequency vector
% interpolate dB on logarithmic frequency axis
H_dB = interp1(log(f_Hz+eps), g_dBc1Hz, log(abs(fb_Hz)+eps), 'linear');
% dB => magnitude
H = 10 .^ (H_dB / 20);
% plot on linear f axis
figure(); hold on;
plot(fftshift(fb_Hz), fftshift(mag2dB(H)), 'k');
xlabel('f/Hz'); ylabel('dBc in 1 Hz');
title('phase noise spectrum (linear frequency axis)');
% plot on logarithmic f axis
figure(); hold on;
ix = find(fb_Hz > 0);
semilogx(fb_Hz(ix), mag2dB(H(ix)), 'k');
xlabel('f/Hz'); ylabel('dBc in 1 Hz');
title('phase noise spectrum (logarithmic frequency axis)');
% apply frequency response to signal
sig = sig .* H;
% add carrier
sig (1) = A;
% convert to time domain
td = ifft(sig);
% determine phase
% for an ideal carrier, it would be zero
% any deviation from zero is phase error
ph_deg = angle(td) * 180 / pi;
figure();
ix = 1:1000;
plot(ix / n / deltaF_Hz, ph_deg(ix));
xlabel('time / s');
ylabel('phase error / degrees');
title('phase of example signal (first 1000 samples)');
figure();
hist(ph_deg, 300);
title('histogram of example signal phase');
xlabel('phase error / degrees');
% RMS average
% note: exp(1i*x) ~ x
% as with magnitude, power/energy is proportional to phase error squared
r = RMS(ph_deg);
% knowing that the signal does not have an amplitude component, we can alternatively
% determine the RMS phase error from the power of the phase noise component
% exclude carrier:
sig(1) = 0;
% 3rd alternative to calculate RMS phase error
rAlt_deg = sqrt(sig * sig') * 180 / pi
end
% gets a vector of frequencies corresponding to n FFT bins, when the
% frequency spacing between two adjacent bins is deltaF_Hz
function fb_Hz = FFT_freqbase(n, deltaF_Hz)
fb_Hz = 0:(n - 1);
fb_Hz = fb_Hz + floor(n / 2);
fb_Hz = mod(fb_Hz, n);
fb_Hz = fb_Hz - floor(n / 2);
fb_Hz = fb_Hz * deltaF_Hz;
end
% Root-mean-square average
function r = RMS(vec)
r = sqrt(vec * vec' / numel(vec));
end
% convert a magnitude to dB
function r = mag2dB(vec)
r = 20*log10(abs(vec) + eps);
end
# Instantiate the SIIR object. Pass the cutoff frequency
# Fc and the sample rate Fs in Hz. Also define the input
# and output fixed-point type. W=(wl, iwl) where
# wl = word-length and iwl = integer word-length. This
# example uses 23 fraction bits and 1 sign bit.
>>> from siir import SIIR
>>> flt = SIIR(Fstop=1333, Fs=48000, W=(24,0))
# Plot the response of the fixed-point coefficients
>>> plot(flt.hz, 20*log10(flt.h)
# Create a testbench and run a simulation
# (get the simulated response)
>>> from myhdl import Simulation
>>> tb = flt.TestFreqResponse(Nloops=128, Nfft=1024)
>>> Simulation(tb).run()
>>> flt.PlotResponse()
# Happy with results generate the Verilog and VHDL
>>> flt.Convert()
% *******************************************************
% 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
//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
