% ****************************************************************
% find sub-sample delay and scaling factor between two cyclic signals
% to maximize crosscorrelation
% Markus Nentwig, 120609_v1.1
% ****************************************************************
function iterDelayEstDemo();
close all;
n = 1024;
% n = 1024 * 256; disp('*** test: long signal enabled ***');
% ****************************************************************
% random signal
% ****************************************************************
fd = randn(1, n) + 1i * randn(1, n);
% ****************************************************************
% lowpass filter
% ****************************************************************
f = binFreq(n);
fd(abs(f) > 0.045) = 0;
s1 = real(ifft(fd)) * sqrt(n);
% ****************************************************************
% create delayed 2nd signal
% ****************************************************************
dTest_samples = 12.3456;
cTest = 1.23456;
% cTest = cTest + 1i; disp('*** test: complex coeff enabled ***');
% cTest = -cTest; disp('*** test: negative coeff enabled ***');
s2 = cTest * cs_delay(s1, 1, dTest_samples);
% s2 = s2 + 0.5*randn(size(s2)); disp('*** test: noise enabled ***');
% ****************************************************************
% estimate delay
% ****************************************************************
[delay_samples, coeff] = iterDelayEst(s1, s2);
% ****************************************************************
% correct it
% ****************************************************************
s2a = cs_delay(s1, 1, delay_samples);
s2b = s2a * coeff;
figure(); hold on;
h = plot(real(s1), 'k'); set(h, 'lineWidth', 3);
h = plot(real(s2), 'b'); set(h, 'lineWidth', 3);
h = plot(real(s2a), 'r'); set(h, 'lineWidth', 1);
h = plot(real(s2b), 'm'); set(h, 'lineWidth', 1);
xlim([1, numel(s1)]);
xlabel('samples');
legend('s1', 's2', 's2 un-delayed', 's2 un-delayed and scaled');
title('test signals');
format long;
disp('nominal delay of s2 relative to s1')';
dTest_samples
disp('iterDelayEst() returned:');
delay_samples
disp('original scaling factor:');
cTest
disp('estimated scaling factor:');
coeff
end
% ****************************************************************
% estimates delay and scaling factor
% ****************************************************************
function [delay_samples, coeff] = iterDelayEst(s1, s2)
s1 = s1(:) .'; % force row vectors
s2 = s2(:) .';
rflag = isreal(s1) && isreal(s2);
n = numel(s1);
halfN = floor(n/2);
assert(numel(s2) == n, 'signals must have same length');
% ****************************************************************
% constants
% ****************************************************************
% exit if uncertainty below threshold
thr_samples = 1e-7;
% exit after fixed number of iterations
nIter = 25;
% frequency domain representation of signals
fd1 = fft(s1);
fd2 = fft(s2);
% first round: No delay was applied
tau = [];
fd2Tau = fd2; % delayed s2 in freq. domain
% frequency corresponding to each FFT bin -0.5..0.5
f = binFreq(n);
% uncertainty plot data
e = [];
% normalization factor
nf = sqrt((fd1 * fd1') * (fd2 * fd2')) / n; % normalizes to 1
% search window:
% known maximum and two surrounding points
x1 = -1;
x2 = -1;
x3 = -1;
y1 = -1;
y2 = -1;
y3 = -1;
% ****************************************************************
% iteration loop
% ****************************************************************
for count = 1:nIter
% ****************************************************************
% crosscorrelation with time-shifted signal
% ****************************************************************
xcorr = abs(ifft(fd2Tau .* conj(fd1)))/ nf;
% ****************************************************************
% detect peak
% ****************************************************************
if isempty(tau)
% ****************************************************************
% startup
% initialize with three adjacent bins around peak
% ****************************************************************
ix = find(xcorr == max(xcorr));
ix = ix(1); % use any, if multiple bitwise identical peaks
% indices of three bins around peak
ixLow = mod(ix-1-1, n) + 1; % one below
ixMid = ix;
ixHigh = mod(ix-1+1, n) + 1; % one above
% delay corresponding to the three bins
tauLow = mod(ixLow -1 + halfN, n) - halfN;
tauMid = mod(ixMid -1 + halfN, n) - halfN;
tauHigh = mod(ixHigh -1 + halfN, n) - halfN;
% crosscorrelation value for the three bins
xcLow = xcorr(ixLow);
xcMid = xcorr(ixMid);
xcHigh = xcorr(ixHigh);
x1 = tauLow;
x2 = tauMid;
x3 = tauHigh;
y1 = xcLow;
y2 = xcMid;
y3 = xcHigh;
else
% ****************************************************************
% only main peak at first bin is of interest
% ****************************************************************
tauMid = tau;
xcMid = xcorr(1);
if xcMid > y2
% ****************************************************************
% improve midpoint
% ****************************************************************
if tauMid > x2
% midpoint becomes lower point
x1 = x2;
y1 = y2;
else
% midpoint becomes upper point
x3 = x2;
y3 = y2;
end
x2 = tauMid;
y2 = xcMid;
elseif tauMid < x2
% ****************************************************************
% improve low point
% ****************************************************************
assert(tauMid >= x1); % bitwise identical is OK
assert(tauMid > x1 || xcMid > y1); % expect improvement
x1 = tauMid;
y1 = xcMid;
elseif tauMid > x2
% ****************************************************************
% improve high point
% ****************************************************************
assert(tauMid <= x3); % bitwise identical is OK
assert((tauMid < x3) || (xcMid > y3)); % expect improvement
x3 = tauMid;
y3 = xcMid;
else
assert(false, '?? evaluated for existing tau ??');
end
end
% ****************************************************************
% calculate uncertainty (window width)
% ****************************************************************
eIter = abs(x3 - x1);
e = [e, eIter];
if eIter < thr_samples
% disp('threshold reached, exiting');
break;
end
if y1 == y2 || y2 == y3
% reached limit of numerical accuracy on one side
usePoly = 0;
else
% ****************************************************************
% fit 2nd order polynomial and find maximum
% ****************************************************************
num = (x2^2-x1^2)*y3+(x1^2-x3^2)*y2+(x3^2-x2^2)*y1;
denom = (2*x2-2*x1)*y3+(2*x1-2*x3)*y2+(2*x3-2*x2)*y1;
if denom ~= 0
tau = num / denom;
% is the point within [x1, x3]?
usePoly = ((tau > x1) && (tau < x3));
else
usePoly = 0;
end
end
if ~usePoly
% revert to linear interpolation on the side with the
% less-accurate outer sample
% Note: There is no guarantee that the side with the more accurate
% outer sample is the right one, as the samples aren't
% placed on a regular grid!
% Therefore, iterate to improve the "worse" side, which will
% eventually become the "better side", and iteration converges.
tauLow = (x1 + x2) / 2;
tauHigh = (x2 + x3) / 2;
if y1 < y3
o = [tauLow, tauHigh];
else
o = [tauHigh, tauLow];
end
% don't choose point that is identical to one that is already known
tau = o(1);
if tau == x1 || tau == x2 || tau == x3
tau = o(2);
if tau == x1 || tau == x2 || tau == x3
break;
end
end
end
% ****************************************************************
% advance 2nd signal according to location of maximum
% phase shift in frequency domain - delay in time domain
% ****************************************************************
fd2Tau = fd2 .* exp(2i * pi * f * tau);
end % for
% ****************************************************************
% plot the uncertainty (window size) over the number of iterations
% ****************************************************************
if true
figure(); semilogy(e, '+-'); grid on;
xlabel('iteration');
title('uncertainty in delay');
end
% ****************************************************************
% the delay estimate is the final location of the delay that
% maximized crosscorrelation (center of window).
% ****************************************************************
delay_samples = x2;
% ****************************************************************
% Coefficient: Turn signal 1 into signal 2
% ****************************************************************
coeff = fd2Tau * fd1' ./ (fd1 * fd1');
% ****************************************************************
% chop roundoff error, if input signals are known to be
% real-valued.
% ****************************************************************
if rflag
coeff = real(coeff);
end
end
% ****************************************************************
% frequency corresponding to FFT bin
% ****************************************************************
function f = binFreq(n)
f = (mod(((0:n-1)+floor(n/2)), n)-floor(n/2))/n;
end
% ****************************************************************
% delay by phase shift
% needed only for demo code
% ****************************************************************
function waveform = cs_delay(waveform, rate_Hz, delay_s)
rflag = isreal(waveform);
n = numel(waveform);
cycLen_s = n / rate_Hz;
nCyc = delay_s / cycLen_s();
f = 0:(n - 1);
f = f + floor(n / 2);
f = mod(f, n);
f = f - floor(n / 2);
phase = -2 * pi * f * nCyc;
rot = exp(1i*phase);
waveform = ifft(fft(waveform) .* rot);
if rflag
waveform = real(waveform);
end
end
NZ = 1; % number of ZEROS in the filter to be designed
NP = 4; % number of POLES in the filter to be designed
NG = 10; % number of gain measurements
fmin = 100; % lowest measurement frequency (Hz)
fmax = 3000; % highest measurement frequency (Hz)
fs = 10000; % discrete-time sampling rate
Nfft = 512; % FFT size to use
df = (fmax/fmin)^(1/(NG-1)); % uniform log-freq spacing
f = fmin * df .^ (0:NG-1); % measurement frequency axis
% Gain measurements (synthetic example = triangular amp response):
Gdb = 10*[1:NG/2,NG/2:-1:1]/(NG/2); % between 0 and 10 dB gain
% Must decide on a dc value.
% Either use what is known to be true or pick something "maximally
% smooth". Here we do a simple linear extrapolation:
dc_amp = Gdb(1) - f(1)*(Gdb(2)-Gdb(1))/(f(2)-f(1));
% Must also decide on a value at half the sampling rate.
% Use either a realistic estimate or something "maximally smooth".
% Here we do a simple linear extrapolation. While zeroing it
% is appealing, we do not want any zeros on the unit circle here.
Gdb_last_slope = (Gdb(NG) - Gdb(NG-1)) / (f(NG) - f(NG-1));
nyq_amp = Gdb(NG) + Gdb_last_slope * (fs/2 - f(NG));
Gdbe = [dc_amp, Gdb, nyq_amp];
fe = [0,f,fs/2];
NGe = NG+2;
% Resample to a uniform frequency grid, as required by ifft.
% We do this by fitting cubic splines evaluated on the fft grid:
Gdbei = spline(fe,Gdbe); % say `help spline'
fk = fs*[0:Nfft/2]/Nfft; % fft frequency grid (nonneg freqs)
Gdbfk = ppval(Gdbei,fk); % Uniformly resampled amp-resp
figure(1);
semilogx(fk(2:end-1),Gdbfk(2:end-1),'-k'); grid('on');
axis([fmin/2 fmax*2 -3 11]);
hold('on'); semilogx(f,Gdb,'ok');
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
title(['Measured and Extrapolated/Interpolated/Resampled ',...
'Amplitude Response']);
Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies
S = 10 .^ (Sdb/20); % convert to linear magnitude
s = ifft(S); % desired impulse response
s = real(s); % any imaginary part is quantization noise
tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
disp(sprintf(['Time-limitedness check: Outer 20%% of impulse ' ...
'response is %0.2f %% of total rms'],tlerr));
% = 0.02 percent
if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
error('Increase Nfft and/or smooth Sdb');
end
figure(2);
plot(s, '-k'); grid('on'); title('Impulse Response');
xlabel('Time (samples)'); ylabel('Amplitude');
c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
% Check aliasing of cepstrum (in theory there is always some):
caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
disp(sprintf(['Cepstral time-aliasing check: Outer 20%% of ' ...
'cepstrum holds %0.2f %% of total rms'],caliaserr));
% = 0.09 percent
if caliaserr>1.0 % arbitrary limit
error('Increase Nfft and/or smooth Sdb to shorten cepstrum');
end
% Fold cepstrum to reflect non-min-phase zeros inside unit circle:
% If complex:
% cf = [c(1), c(2:Ns-1)+conj(c(Nfft:-1:Ns+1)), c(Ns), zeros(1,Nfft-Ns)];
cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
Cf = fft(cf); % = dB_magnitude + j * minimum_phase
Smp = 10 .^ (Cf/20); % minimum-phase spectrum
Smpp = Smp(1:Ns); % nonnegative-frequency portion
wt = 1 ./ (fk+1); % typical weight fn for audio
wk = 2*pi*fk/fs;
[B,A] = invfreqz(Smpp,wk,NZ,NP,wt);
Hh = freqz(B,A,Ns);
figure(3);
plot(fk,db([Smpp(:),Hh(:)])); grid('on');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Frequency Response');
legend('Desired','Filter');
function [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB)
%
% SNR_Set(x, Desired_SNR_dB) returns the real-valued
% input 'Signal' contaminated with normally-distributed,
% zero-mean, random noise. The signal-to-noise ratio
% (SNR in dB) of the output 'Noisy_Signal' signal is
% controlled by the input 'Desired_SNR_dB' argument measured
% in dB.
% Example:
% Npts = 128; % Number of time samples
% n = 0:Npts-1; % Time-domain index
% Signal = 3*sin(2*pi*3*n/Npts); % Real-valued signal
% Desired_SNR_dB = 3; % Set SNR of output 'Noisy_Signal' to +3 dB
% [Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB);
%
% Author: Richard Lyons [December 2011]
%******************************************
Npts = length(Signal); % Number of input time samples
Noise = randn(1,Npts); % Generate initial noise; mean zero, variance one
Signal_Power = sum(abs(Signal).*abs(Signal))/Npts;
Noise_Power = sum(abs(Noise).*abs(Noise))/Npts;
%Initial_SNR = 10*(log10(Signal_Power/Noise_Power));
K = (Signal_Power/Noise_Power)*10^(-Desired_SNR_dB/10); % Scale factor
New_Noise = sqrt(K)*Noise; % Change Noise level
%New_Noise_Power = sum(abs(New_Noise).*abs(New_Noise))/Npts
%New_SNR = 10*(log10(Signal_Power/New_Noise_Power))
Noisy_Signal = Signal + New_Noise;
'SNR_Set()' Function Test Code:
% Filename SNR_Set_test.m
%
% Tests the 'SNR_Set()" function. Adds a predefined
% amount of random noise to a noise-free signal such that
% the noisy signal has a desired signal-to-noise ratio (SNR).
%
% Author: Richard Lyons [December 2011]
clear, clc
% Create a noise-free signal
Npts = 128; % Number of time samples
n = 0:Npts-1; % Time-domain index
Cycles = 5; % Integer number of cycles in noise-free sinwave signal
Signal = 3*sin(2*pi*Cycles*n/Npts); % Real-valued noise-free signal
Desired_SNR_dB = 3 % Set desired SNR in dB
[Noisy_Signal] = SNR_Set(Signal, Desired_SNR_dB); % Generate noisy signal
% Plot original and 'noisy' signals
figure(1)
subplot(2,1,1)
plot(n, Signal, '-bo', 'markersize', 2)
title('Original Signal')
grid on, zoom on
subplot(2,1,2)
plot(n, Noisy_Signal, '-bo', 'markersize', 2)
title('Noisy Signal')
xlabel('Time-samples')
grid on, zoom on
% Measure SNR in freq domain
Spec = fft(Noisy_Signal);
Spec_Mag = abs(Spec); % Spectral magnitude
figure(2)
plot(Spec_Mag, '-bo', 'markersize', 2)
title('Spec Mag of Noisy Signal')
xlabel('Freq-samples'), ylabel('Linear')
grid on, zoom on
Signal_Power = Spec_Mag(Cycles+1)^2 + Spec_Mag(Npts-Cycles+1)^2;
Noise_Power = sum(Spec_Mag.^2) -Signal_Power;
Measured_SNR = 10*log10(Signal_Power/Noise_Power)
#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
#define pi 3.14159265
// class complex
class complex
{
public:
//This class member function set the values of the class members real and imag
void setvals(double r, double i)
{
real=r;
imag=i;
}
//This function get the values (by reference) of the class members real and imag
void getvals(double &r, double &i)
{
r=real;
i=imag;
}
private:
/* The class members real and imag represents the real and imaginary parts
of a complex number, respectively.*/
double real;
double imag;
}; // end of class complex
//**********************************************************************************************
/* function dft receives two parameters x and N. x is the input sequence containing
the discrete-time samples and N is the number of DFT points. The complex type paramater
xk receives a sequence of N-points frequency samples (this is the DFT of the input sequence*/
void dft(double x[], int N, complex xk[])
{
double r,i;
for(int k=0;k<N;k++) //k is the index of the N-point frequency sequence
{
xk[k].setvals(0,0); //initialize the values of the counters to zero
for(int n=0;n<N;n++) //n is the index of the discrete-time sequence
{
xk[k].getvals(r,i); // get the values of real and imag
/*This is the computation of the real and imaginary parts of the DFT.
The class member function setvals is used to update the value of the accumulators
contaning the real and imaginary parts of the DFT samples*/
xk[k].setvals(r + x[n]*cos(2*pi*k/N*n),i + x[n]*sin(2*pi*k/N*n));
} //end inner for-loop
} //end outer for-loop
}//end dft function
//************************************************************************************************
//Declaration of the main() function
int main()
{
int size; // The number of DFT samples
double r,i; // r = real part of the DFT; i = imaginary part of the DFT
cout<<"Enter the desired number of DFT samples: ";
cin>>size;
double *x = new double[size]; //dynamic array x represent the discrete-time samples
complex *xk = new complex[size]; // dynamic array xk represent the DFT samples
cout<<"\nEnter the sequence x[n]: ";
for(int h=0;h<size;h++)
{
cin>>x[h];
}
//Computation of the DFT of x
dft(x, size, xk);
system("CLS"); //clear screen
cout<<"Discrete Fourier Transform: "<<endl<<endl;
for(int h=0;h<size;h++)
{
cout.precision(3); // display three (3) decimal places
xk[h].getvals(r,i); //get the DFT samples values
cout<<"xk["<<h<<"] = "<<fixed<< r <<" + "<<fixed<< i <<"j"<<endl; // print the DFT samples
}
cout<<endl;
return 0;
}// end main function
% The function NDFT2 computes the sinc pattern of a linear array of
% sensors. The function receives three input: the interelement spacing d,
% the array weigths a, and the zero padding desired Npad.
function NDFT2(d, a, Npad)
k = -Npad/2 : Npad/2-1; % index
u = pi*k/Npad; % u is a vector with elements in the range of -pi/2 to pi/2
uk = sin(u);
n = 0:Npad-1; % n is an index
m = 1:Npad; % m is an index
% Wavenumber K = 2*pi/landa
% d = is a fraction or a multiple of lambda.
% v = K*d*uk(m).
v = 2*pi*d*uk(m)';
Dn(m,n+1) = exp(j*v*n);
puk = Dn * a'; % Computing the Beampattern
% Plot of the Beampatterns
figure(1); subplot(2,1,1); plot(uk,20*log10(abs(puk)));grid on; xlabel('sin(theta)'); ylabel('Magnitude (dB)')
axis([-1 1 -100 0])
subplot(2,1,2); plot(uk,puk);grid on; xlabel('sin(theta)'); ylabel('Magnitude')
axis([-1 1 -1 1])
warning off;
% Plot the beampattern as a polar graph
figure(2); polar(u',puk); hold on; polar(-u',puk);hold off
% Function pattern()
%
% The function pattern() computes and plots the beampattern of a
% Linear Array of sensors. This function has three inputs: the number of elements in
% the array, the pointing direction of the beam pattern (an angular
% direction in radians), and a constant spacing between the elements of
% the array (as fraction of the wavelenght(lambda)of the transmitted
% signal. The optimum spacing between elements of the array is 0.5.
% You must also select any of the windows presented in the menu.
% Windowing techniques are used to reduced the sidelobes of the pattern.
% The list of available windows is the following:
%
% @bartlett - Bartlett window.
% @barthannwin - Modified Bartlett-Hanning window.
% @blackman - Blackman window.
% @blackmanharris - Minimum 4-term Blackman-Harris window.
% @bohmanwin - Bohman window.
% @chebwin - Chebyshev window.
% @flattopwin - Flat Top window.
% @gausswin - Gaussian window.
% @hamming - Hamming window.
% @hann - Hann window.
% @kaiser - Kaiser window.
% @nuttallwin - Nuttall defined minimum 4-term Blackman-Harris window.
% @parzenwin - Parzen (de la Valle-Poussin) window.
% @rectwin - Rectangular window.
% @tukeywin - Tukey window.
% @triang - Triangular window.
%
% For example: pattern(21, pi/4, 0.5);
% Plots the beampattern of an linear array with 21 elements equally spaced
% at a half of the wavelenght(lambda/2), and a pointing
% direction of 45 degrees. For uniform arrays use a rectangular
% window (rectwin).
function pattern(array_number, angular_direction, spacing)
close all
clc
N = array_number;
delta = angular_direction;
d = spacing;
Npad=1024;
n=0:N-1;
delta = 2*pi*d*sin(delta);
an=1/N*exp(j*n*delta);
for i=0:500000,
option = menu('Choose the desired Window', 'Bartlett', 'Barthannwin', 'Blackman', 'Blackmanharris', 'Bohmanwin', 'Chebwin', 'Flattopwin', 'Gausswin', 'Hamming', 'Hann', 'Kaiser', 'Nuttallwin', 'Parzenwin', 'Rectwin', 'Tukeywin', 'Triang', 'Exit');
switch option
case 1
close all
clear a;
a=an;
a = a.*bartlett(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 2
close all
clear a;
a=an;
a = a.*barthannwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 3
close all
clear a;
a=an;
a = a.*blackman(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 4
close all
clear a;
a=an;
a = a.*blackmanharris(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 5
close all
clear a;
a=an;
a = a.*bohmanwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 6
close all
clear a;
a=an;
a = a.*chebwin(N,40)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 7
close all
clear a;
a=an;
a = a.*flattopwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 8
close all
clear a;
a=an;
a = a.*gausswin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 9
close all
clear a;
a=an;
a = a.*hamming(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 10
close all
clear a;
a=an;
a = a.*hann(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 11
close all
clear a;
a=an;
a = a.*kaiser(N,1)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 12
close all
clear a;
a=an;
a = a.*nuttallwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 13
close all
clear a;
a=an;
a = a.*parzenwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 14
close all
clear a;
a=an;
a = a.*rectwin(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 15
close all
clear a;
a=an;
a = a.*tukeywin(N,0)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 16
close all
clear a;
a=an;
a = a.*triang(N)';
a(Npad)=0;
NDFT2(d, a, Npad);
case 17
break;
end
end
% function echo_filter:
% 1) computes and returns the impulse response h of the LTI system (echo filter)
% 2) Filter the input signal. The output signal is yecho.
%
% INPUTS:
% signal = input signal
% times = vector containing the time delays (in seconds) of the repetitions(echoes)
% attenuations = vector containing the respective attenuations of the echoes
% fs = sampling frequency of the input signal.
%
% OUTPUTS:
% yecho = output signal (passed through the echo filter)
% h = impulse response of the echo filter.
function [yecho, h] = echo_filter(signal, times, attenuations, fs)
h(1)=1; % This is the first coefficient of the echo filter
% Computing the impulse response h
for i=1:length(times),
samples(i) = times(i)*fs; % Calculating the sample-times (N = t*fs)
h(floor(samples(i))) = attenuations(i); % Impulse response coeficients
end
% #########################################################################
% You may use the following implementation instead of the illustrated
% above:
% samples = times*fs;
% h(floor(samples)) = attenuations;
% In this case, the implementation is done without a for-loop.
% #########################################################################
yecho = filter(h,1,signal(:,1)); % filtering of the input signal results in a signal with echoes
% You may test the filter using the Matlab signal "splat" (write load splat
% in the Matlab command window). Use function sound() in order to
% listening both, the input and the output, signals.
% ******************************************************
% Symbol error rate of Quadrature Amplitude Modulation
% ******************************************************
% uncoded, additive white Gaussian noise channel
% Implements [1] Figure 5.2-16 on a wider SNR range
% [1] John G. Proakis
% "Digital Communications"
% 4th edition, McGraw-Hill", 2001
% http://www.amazon.com/Digital-Communications-John-Proakis/dp/0072321113
% Note that different definitions are in common use:
% a) {symbol, bit} errors
% b) signal-to-noise ratio per {symbol, bit}
close all; clear all;
SNR_perSymbol_dB=[-3:0.1:42];
% definition: Below [1] 5.2-78
SNR_perSymbol=10.^(SNR_perSymbol_dB/10);
plotPerBitX=[]; plotY=[];
% number of constellation points (M-ary QAM)
for M=[4, 16, 64, 256, 1024]
% the energy of each symbol is used to transmit log2(M) bits
SNR_perBit=SNR_perSymbol/log2(M);
% [1] 5.2-78 argument of Q(...)
Qarg=sqrt(3/(M-1)*SNR_perSymbol);
% [1] 2.1-98
Q=1/2*erfc(Qarg/sqrt(2));
% [1] eq. 5.2-77
% probability of error for PAM per quadrature signal
PsqrtM=2*(1-1/sqrt(M))*Q;
% [1] eq. 5.2-79
PM=1-(1-PsqrtM).^2;
plotPerBitX=[plotPerBitX; 10*log10(SNR_perBit)];
plotY=[plotY; PM];
end
figure(1);
h = semilogy(SNR_perSymbol_dB', plotY'); grid on;
set(h, 'lineWidth', 3);
title('M-ary QAM'); xlabel('SNR per symbol (dB)'); ylabel('symbol error rate');
xlim([0, 40]); ylim([1e-6, 1]);
legend({'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'});
figure(2);
h = semilogy(plotPerBitX', plotY'); grid on; ylabel('symbol error rate');
set(h, 'lineWidth', 3);
title('M-ary QAM'); xlabel('SNR per bit (dB)');
xlim([-3, 30]); ylim([1e-6, 1]);
legend({'QPSK', '16QAM', '64QAM', '256QAM', '1024QAM'});
%Function to calculate the Average of a Voice signal
function [y] = pow_1(x)
N = length(x); % Length of voice signal 'x'
xold =0.0; %initialize it to zero
for n = 1:N
sumx = xold+(x(n)^2);
xold = sumx;
end
y = sumx/N
// To Design an Analog Butterworth Filter
//For the given cutoff frequency and filter order
//Wc = 500 Hz
omegap = 500; //pass band edge frequency
omegas = 1000;//stop band edge frequency
delta1_in_dB = -3;//PassBand Ripple in dB
delta2_in_dB = -40;//StopBand Ripple in dB
delta1 = 10^(delta1_in_dB/20)
delta2 = 10^(delta2_in_dB/20)
//Caculation of filter order
N = log10((1/(delta2^2))-1)/(2*log10(omegas/omegap))
N = ceil(N) //Rounding off nearest integer
omegac = omegap;
h=buttmag(N,omegac,1:1000);//Analog Butterworth filter magnitude response
mag=20*log10(h);//Magntitude Response in dB
plot2d((1:1000),mag,[0,-180,1000,20]);
a=gca();
a.thickness = 3;
a.foreground = 1;
a.font_style = 9;
xgrid(5)
xtitle('Magnitude Response of Butterworth LPF Filter Cutoff frequency = 500 Hz','Analog frequency in Hz--->','Magnitude in dB -->');
//Design of FIR Filter using Frquency Sampling Technique
//Low Pass Filter Design
//Cutoff Frequency Wc = pi/2
//M = Filter Lenth = 7
clear;
clc;
M = 7;
N = ((M-1)/2)+1
wc = %pi/2;
for k =1:M
w(k) = ((2*%pi)/M)*(k-1);
if (w(k)>=wc)
k-1
break
end
end
for i = 1:k-1
Hr(i) = 1;
G(i) = ((-1)^(i-1))*Hr(i);
end
for i = k:N
Hr(i) = 0;
G(i) = ((-1)^(i-1))*Hr(i);
end
h = zeros(1,M);
for n = 1:M
for k = 2:N
h(n) = G(k)*cos((2*%pi/M)*(k-1)*((n-1)+(1/2)))+h(n);
end
h(n) = (1/M)* (G(1)+2*h(n));
end
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xtitle('Frequency Response of LPF with Normalized cutoff =0.5','Normalized Frequency W ------>','Magnitude Response H(w)---->');
xgrid(1);
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//Band Stop FILTER (or)Band Reject Filter
clear all;
clc;
close;
M = 7 //Filter length = 11
Wc = [%pi/5,3*%pi/4]; //Digital Cutoff frequency
Wc2 = Wc(2)
Wc1 = Wc(1)
Tuo = (M-1)/2 //Center Value
hd = zeros(1,M);
W = zeros(1,M);
for n = 1:M
if (n == Tuo+1)
hd(n) = 1-((Wc2-Wc1)/%pi);
else hd(n)=(sin(%pi*((n-1)-Tuo))-sin(Wc2*((n-1)-Tuo))+sin(Wc1*((n-1)-Tuo)))/(((n-1)-Tuo)*%pi);
end
if(abs(hd(n))<(0.00001))
hd(n)=0;
end
end
hd
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp(h,'Filter Coefficients are')
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR BPF using Rectangular window M=7')
xgrid(1)
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//Band PASS FILTER - Window Based
clear all;
clc;
close;
M = 7 //Filter length = 7
Wc = [%pi/5,3*%pi/4]; //Digital Cutoff frequency
Wc2 = Wc(2)
Wc1 = Wc(1)
Tuo = (M-1)/2 //Center Value
hd = zeros(1,M);
W = zeros(1,M);
for n = 1:M
if (n == Tuo+1)
hd(n) = (Wc2-Wc1)/%pi;
else
n
hd(n) = (sin(Wc2*((n-1)-Tuo)) -sin(Wc1*((n-1)-Tuo)))/(((n-1)-Tuo)*%pi);
end
if(abs(hd(n))<(0.00001))
hd(n)=0;
end
end
hd;
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp(h,'Filter Coefficients are')
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(2*fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR BPF using Rectangular window M=7')
xgrid(1)
//PROGRAM TO DESIGN AND OBTAIN THE FREQUENCY RESPONSE OF FIR FILTER
//HIGH PASS FILTER -WINDOW BASED
clear all;
clc;
close;
M = 7 //Filter length = 7
Wc = %pi/4; //Digital Cutoff frequency
Tuo = (M-1)/2 //Center Value
for n = 1:M
if (n == Tuo+1)
hd(n) = 1-Wc/%pi;
else
hd(n) = (sin(%pi*((n-1)-Tuo)) -sin(Wc*((n-1)-Tuo)))/(((n-1)-Tuo)*%pi);
end
end
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp('Filter Coefficients are')
h;
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR HPF using Rectangular window M=7')
//Program to design a FIR Low Pass Filter- Window Based
//Technique
clear all;
clc;
close;
M = 7 //Filter length = 7
Wc = %pi/4; //Digital Cutoff frequency
Tuo = (M-1)/2 //Center Value
for n = 1:M
if (n == Tuo+1)
hd(n) = Wc/%pi;
else
hd(n) = sin(Wc*((n-1)-Tuo))/(((n-1)-Tuo)*%pi);
end
end
//Rectangular Window
for n = 1:M
W(n) = 1;
end
//Windowing Fitler Coefficients
h = hd.*W;
disp('Filter Coefficients are')
h;
[hzm,fr]=frmag(h,256);
hzm_dB = 20*log10(hzm)./max(hzm);
plot(fr,hzm_dB)
xlabel('Normalized Digital Frequency W');
ylabel('Magnitude in dB');
title('Frequency Response 0f FIR LPF using Rectangular window M=7')
//Band Pass FIlter of length M = 16
//Lower Cutoff frequency fp = 0.2 and Upper Cutoff frequency fs = 0.3
// Choose the number of cosine functions and create a dense grid
// in [0,0.1) and [0.2,0.35] and [0.425,0.5]
//magnitude for pass band = 1 & stop band = 0 (i.e) [0 1 0]
//Weighting function =[10 1 10]
clear all;
clc;
close;
hn = 0;
hm = 0;
hn=eqfir(16,[0 .1;.2 .35;.425 .5],[0 1 0],[10 1 10]);
[hm,fr]=frmag(hn,256);
disp(hn,'The Filter Coefficients are:')
figure
plot(.5*(0:255)/256,20*log10(frmag(hn,256)));
a = gca();
xlabel('Normalized Digital Frequency fr');
ylabel('Magnitude in dB');
title('Frequency Response of FIR BPF using REMEZ algorithm M=16')
xgrid(2)
//Low Pass FIlter of length M = 11
//Pass band Edge frequency fp = 0.1 and a Stop edge frequency fs = 0.15
// Choose the number of cosine functions and create a dense grid
// in [0,0.2) and [0.25,0.5)
//magnitude for pass band = 1 & stop band = 0 (i.e) [1 0]
//Weighting function =[1 1]
clear all;
clc;
close;
M =11;
hn=eqfir(11,[0,0.2;0.25,0.5],[1 0],[1 1]);
[hm,fr]=frmag(hn,256);
disp(hn,'The Filter Coefficients are:')
figure
plot(.5*(0:255)/256,20*log10(frmag(hn,256)));
xlabel('Normalized Digital Frequency fr');
ylabel('Magnitude in dB');
title('Frequency Response of FIR LPF using REMEZ algorithm M=11')
xgrid(2)
//Program To convert analog IIR filter into digital IIR filter using
//Bilinear Transformation
clear all;
clc;
close;
s = poly(0,'s');
H = (s+0.1)/((s+0.1)^2+16);
T = 0.5;
z = poly(0,'z');
Hz = horner(H,(2/T)*((z-1)/(z+1)))
