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
% Filename: FFT_Twiddles_Find_DSPrelated.m
% Computes 'Decimation in Frequency' or 'Decimation
% in Time' Butterfly twiddle factors, for radix-2 FFTs
% with in-order input indices and scrambled output indices.
%
% To use, do two things: (1) define FFT size 'N'; and
% (2) define the desired 'Structure', near line 17-18,
% as 'Dec_in_Time' or 'Dec_in_Freq'.
%
% Author: Richard Lyons, November, 2011
%******************************************
clear, clc
% Define input parameters
N = 8; % FFT size (Must be an integer power of 2)
Structure = 'Dec_in_Time'; % Choose Dec-in-time butterflies
%Structure = 'Dec_in_Freq'; % Choose Dec-in-frequency butterflies
% Start of processing
Num_Stages = log2(N); % Number of stages
StageStart = 1; % First stage to compute
StageStop = Num_Stages; % Last stage to compute
ButterStart = 1; %First butterfly to compute
ButterStop = N/2; %Last butterfly to compute
Pointer = 0; %Init 'results' row pointer
for Stage_Num = StageStart:StageStop
if Structure == 'Dec_in_Time'
for Butter_Num = ButterStart:ButterStop
Twid = floor((2^Stage_Num*(Butter_Num-1))/N);
% Compute bit reversal of Twid
Twid_Bit_Rev = 0;
for I = Num_Stages-2:-1:0
if Twid>=2^I
Twid_Bit_Rev = Twid_Bit_Rev + 2^(Num_Stages-I-2);
Twid = Twid -2^I;
else, end
end %End bit reversal 'I' loop
A1 = Twid_Bit_Rev; %Angle A1
A2 = Twid_Bit_Rev + N/2; %Angle A2
Pointer = Pointer +1;
Results(Pointer,:) = [Stage_Num,Butter_Num,A1,A2];
end
else
for Twiddle_Num = 1:N/2^Stage_Num
Twid = (2^Stage_Num*(Twiddle_Num-1))/2; %Compute integer
Pointer = Pointer +1;
Results(Pointer,:) = [Stage_Num,Twiddle_Num,Twid];
end
end % End 'if'
end % End Stage_Num loop
Results(:,1:3), disp(' Stage# Twid# A'), disp(' ')
% WCDMA channelization codes
% source:
% 3GPP TS 25.213 V10.0.0 section 4.3.1.1
%
% parameters:
% sf: spreading factor
% k: code number
%
% The returned code is a column vector with length sf
%
% this code has not been tested extensively. Please verify
% independently that it does what it promises.
function code=UTRA_FDD_chanCode(sf, k)
persistent flag;
persistent codes;
% * ********************************************
% * Build codebook
% * ********************************************
if isempty(flag)
codes={};
flag=1;
% start of recursion: Identity code for sf=1
codes{1, 1}=1;
for sfi=1:8
sfg=2 ^ sfi;
for kgDest=0:2:sfg-2
kgSrc=kgDest/2;
prev=codes{sfg/2, kgSrc+1};
% generation method:
% - duplicate a lower-sf code
% - duplicate and change sign
codes{sfg, kgDest+1}=[prev prev];
codes{sfg, kgDest+2}=[prev -prev];
end
end
end
% * ********************************************
% * look up the requested code from codebook
% * ********************************************
code=transpose(codes{sf, k+1});
% ################## CUT HERE #########################
% Example use (put this into a separate file)
sf=128; codenum=3;
chanCode=UTRA_FDD_chanCode(sf, codenum);
sig=[1 0 0 1 0 0 ]; % signal, row vector
len=size(sig, 2),
% convolve:
s=chanCode * sig; % now matrix, one row per repetition
len=len*sf;
s=reshape(s, 1, len);
% plot
stem(s);
% **************************************************************
% Vectorized Farrow resampler
% M. Nentwig, 2011
% Note: Uses cyclic signals (wraps around)
% **************************************************************
close all; clear all;
% inData contains input to the resampling process (instead of function arguments)
inData = struct();
% **************************************************************
% example coefficients.
% Each column [c0; c1; c2; c3] describes a polynomial for one tap coefficent in fractional time ft [0, 1]:
% tapCoeff = c0 + c1 * ft + c2 * ft ^ 2 + c3 * ft ^ 3
% Each column corresponds to one tap.
% the matrix size may be changed arbitrarily.
%
% The example filter is based on a 6th order Chebyshev Laplace-domain prototype.
% **************************************************************
inData.cMatrix =[ -8.57738278e-3 7.82989032e-1 7.19303539e+000 6.90955718e+000 -2.62377450e+000 -6.85327127e-1 1.44681608e+000 -8.79147907e-1 7.82633997e-2 1.91318985e-1 -1.88573400e-1 6.91790782e-2 3.07723786e-3 -6.74800912e-3
2.32448021e-1 2.52624309e+000 7.67543936e+000 -8.83951796e+000 -5.49838636e+000 6.07298348e+000 -2.16053205e+000 -7.59142947e-1 1.41269409e+000 -8.17735712e-1 1.98119464e-1 9.15904145e-2 -9.18092030e-2 2.74136108e-2
-1.14183319e+000 6.86126458e+000 -6.86015957e+000 -6.35135894e+000 1.10745051e+001 -3.34847578e+000 -2.22405694e+000 3.14374725e+000 -1.68249886e+000 2.54083065e-1 3.22275037e-1 -3.04794927e-1 1.29393976e-1 -3.32026332e-2
1.67363115e+000 -2.93090391e+000 -1.13549165e+000 5.65274939e+000 -3.60291782e+000 -6.20715544e-1 2.06619782e+000 -1.42159644e+000 3.75075865e-1 1.88433333e-1 -2.64135123e-1 1.47117661e-1 -4.71871047e-2 1.24921920e-2
] / 12.28;
% **************************************************************
% Create example signal
% **************************************************************
nIn = 50; % used for test signal generation only
if true
% complex signal
inData.signal = cos(2*pi*(0:(nIn-1)) / nIn) + cos(2*2*pi*(0:(nIn-1)) / nIn) + 1.5;
else
% impulse response
inData.signal = zeros(1, nIn); inData.signal(1) = 1; %
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 1) = 1; % enable to show constant c in first tap
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(2, 1) = 1; % enable to show linear c in first tap
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(3, 1) = 1; % enable to show quadratic c in first tap
% inData.cMatrix = 0 * inData.cMatrix; inData.cMatrix(1, 2) = 1; % enable to show constant c in second tap
end
% **************************************************************
% Resample to the following number of output samples
% must be integer, otherwise arbitrary
% **************************************************************
inData.nSamplesOut = floor(nIn * 6.28);
% **************************************************************
% Set up Farrow resampling
% **************************************************************
nSamplesIn = size(inData.signal, 2);
nSamplesOut = inData.nSamplesOut;
order = size(inData.cMatrix, 1) - 1; % polynomial order
nTaps = size(inData.cMatrix, 2); % number of input samples that contribute to one output sample (FIR size)
% pointer to the position in the input stream for each output sample (row vector, real numbers), starting at 0
inputIndex = (0:nSamplesOut-1) / nSamplesOut * nSamplesIn;
% split into integer part (0..nSamplesIn - 1) ...
inputIndexIntegerPart = floor(inputIndex);
% ... and fractional part [0, 1[
inputIndexFractionalPart = inputIndex - inputIndexIntegerPart;
% **************************************************************
% Calculate output stream
% First constant term (conventional FIR), then linear, quadratic, cubic, ...
% **************************************************************
outStream = zeros(1, inData.nSamplesOut);
for ixOrder = 0 : order
x = inputIndexFractionalPart .^ ixOrder;
% **************************************************************
% Add the contribution of each tap one-by-one
% **************************************************************
for ixTap = 0 : nTaps - 1
c = inData.cMatrix(ixOrder+1, ixTap+1);
% index of input sample that contributes to output via the current tap
% higher tap index => longer delay => older input sample => smaller data index
dataIx = inputIndexIntegerPart - ixTap;
% wrap around
dataIx = mod(dataIx, nSamplesIn);
% array indexing starts at 1
dataIx = dataIx + 1;
delayed = inData.signal(dataIx);
% for each individual output sample (index in row vector),
% - evaluate f = c(order, tapindex) * fracPart .^ order
% - scale the delayed input sample with f
% - accumulate the contribution of all taps
% this implementation performs the operation for all output samples in parallel (row vector)
outStream = outStream + c * delayed .* x;
end % for ixTap
end % for ixOrder
% **************************************************************
% plot
% **************************************************************
xIn = linspace(0, 1, nSamplesIn + 1); xIn = xIn(1:end-1);
xOut = linspace(0, 1, nSamplesOut + 1); xOut = xOut(1:end-1);
figure(); grid on; hold on;
stem(xIn, inData.signal, 'k+-');
plot(xOut, outStream, 'b+-');
legend('input', 'output');
title('Farrow resampling. Signals are cyclic');
function [Windowed_Spec] = Wind_Flattop(Spec)
% Given an input spectral sequence 'Spec', that is the
% FFT of some time sequence 'x', Wind_Flattop(Spec)
% returns a spectral sequence that is equivalent
% to the FFT of a flat-top windowed version of time
% sequence 'x'. The peak magnitude values of output
% sequence 'Windowed_Spec' can be used to accurately
% estimate the peak amplitudes of sinusoidal components
% in time sequence 'x'.
% Input: 'Spec' (a sequence of complex FFT sample values)
%
% Output: 'Windowed_Spec' (a sequence of complex flat-top
% windowed FFT sample values)
%
% Based on Lyons': "Reducing FFT Scalloping Loss Errors
% Without Multiplication", IEEE Signal Processing Magazine,
% DSP Tips & Tricks column, March, 2011, pp. 112-116.
%
% Richard Lyons [December, 2011]
N = length(Spec);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform freq-domain convolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g_Coeffs = [1, -0.94247, 0.44247];
% Compute first two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(1) = g_Coeffs(3)*Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + Spec(1) ...
+g_Coeffs(2)*Spec(2) + g_Coeffs(3)*Spec(3);
Windowed_Spec(2) = g_Coeffs(3)*Spec(N) ...
+g_Coeffs(2)*Spec(1) + Spec(2) ...
+g_Coeffs(2)*Spec(3) + g_Coeffs(3)*Spec(4);
% Compute last two convolved spec samples using spectral 'wrap-around'
Windowed_Spec(N-1) = g_Coeffs(3)*Spec(N-3) ...
+g_Coeffs(2)*Spec(N-2) + Spec(N-1) ...
+g_Coeffs(2)*Spec(N) + g_Coeffs(3)*Spec(1);
Windowed_Spec(N) = g_Coeffs(3)*Spec(N-2) ...
+g_Coeffs(2)*Spec(N-1) + Spec(N) ...
+g_Coeffs(2)*Spec(1) + g_Coeffs(3)*Spec(2);
% Compute convolved spec samples for the middle of the spectrum
for K = 3:N-2
Windowed_Spec(K) = g_Coeffs(3)*Spec(K-2) ...
+g_Coeffs(2)*Spec(K-1) + Spec(K) ...
+g_Coeffs(2)*Spec(K+1) + g_Coeffs(3)*Spec(K+2);
end % % End of 'Wind_Flattop(Spec)' function
% Lagrange interpolation for resampling
% References:
% [1] A digital signal processing approach to Interpolation
% Ronald W. Schafer and Lawrence R. Rabiner
% Proc. IEEE vol 61, pp.692-702, June 1973
% [2] https://ccrma.stanford.edu/~jos/Interpolation/Lagrange_Interpolation.html
% [3] Splitting the Unit delay: Tools for fractional delay filter design
% T. I. Laakso, V. Valimaki, M. Karjalainen, and U. K. Laine
% IEEE Signal Processing Magazine, vol. 13, no. 1, pp. 30..60, January 1996
function lagrangeResamplingDemo()
originDefinition = 0; % see comment for LagrangeBasisPolynomial() below
% Regarding order: From [1]
% "Indeed, it is easy to show that whenever Q is odd, none of the
% impulse responses corresponding to Lagrange interpolation can have
% linear phase"
% Here, order = Q+1 => odd orders are preferable
order = 3;
% *******************************************************************
% Set up signals
% *******************************************************************
nIn = order + 1;
nOut = nIn * 5 * 20;
% Reference data: from [1] fig. 8, linear-phase type
ref = [-0.032, -0.056, -0.064, -0.048, 0, ...
0.216, 0.448, 0.672, 0.864, 1, ...
0.864, 0.672, 0.448, 0.216, 0, ...
-0.048, -0.064, -0.056, -0.032, 0];
tRef = (1:size(ref, 2)) / size(ref, 2);
% impulse input to obtain impulse response
inData = zeros(1, nIn);
inData(1) = 1;
outData = zeros(1, nOut);
outTime = 0:(nOut-1);
outTimeAtInput = outTime / nOut * nIn;
outTimeAtInputInteger = floor(outTimeAtInput);
outTimeAtInputFractional = outTimeAtInput - outTimeAtInputInteger;
evalFracTime = outTimeAtInputFractional - 0.5 + originDefinition;
% odd-order modification
if mod(order, 2) == 0
% Continuity of the impulse response is achieved when support points are located at
% the intersections between adjacent segments "at +/- 0.5"
% For an even order polynomial (odd number of points), this is only possible with
% an asymmetric impulse response
offset = 0.5;
%offset = -0.5; % alternatively, its mirror image
else
offset = 0;
end
% *******************************************************************
% Apply Lagrange interpolation to input data
% *******************************************************************
for ixTap = 0:order
% ixInSample is the input sample that appears at FIR tap 'ixTap' to contribute
% to the output sample
% Row vector, for all output samples in parallel
ixInSample = mod(outTimeAtInputInteger + ixTap - order, nIn) + 1;
% the value of said input sample, for all output samples in parallel
d = inData(ixInSample);
% Get Lagrange polynomial coefficients of this tap
c = lagrangeBasisPolynomial(order, ixTap, originDefinition + offset);
% Evaluate the Lagrange polynomial, resulting in the time-varying FIR tap weight
cTap = polyval(c, evalFracTime);
% FIR operation: multiply FIR tap weight with input sample and add to
% output sample (all outputs in parallel)
outData = outData + d .* cTap;
end
% *******************************************************************
% Plot
% *******************************************************************
figure(); hold on;
h = plot((0:nOut-1) / nOut, outData, 'b-'); set(h, 'lineWidth', 3);
stem(tRef, ref, 'r'); set(h, 'lineWidth', 3);
legend('impulse response', 'reference result');
title('Lagrange interpolation, impulse response');
end
% Returns the coefficients of a Lagrange basis polynomial
% 1 <= order: Polynomial order
% 0 <= evalIx <= order: index of basis function.
%
% At the set of support points, the basis polynomials evaluate as follows:
% evalIx = 1: [1 0 0 ...]
% evalIx = 2: [0 1 0 ...]
% evalIx = 3: [0 0 1 ...]
%
% The support point are equally spaced.
% Example, using originDefinition=0:
% order = 1: [-0.5 0.5]
% order = 2: [-1 0 1]
% order = 3: [-1.5 -0.5 0.5 1.5]
%
% The area around the mid-point corresponds to -0.5 <= x <= 0.5.
% If a resampler implementation uses by convention 0 <= x <= 1 instead, set
% originDefinition=0.5 to offset
% the polynomial.
function polyCoeff = lagrangeBasisPolynomial(order, evalIx, originDefinition)
assert(evalIx >= 0);
assert(evalIx <= order);
tapLocations = -0.5*(order) + (0:order) + originDefinition;
polyCoeff = [1];
for loopIx = 0:order
if loopIx ~= evalIx
% numerator: places a zero in the polynomial via (x-xTap(k)), with k != evalIx
% denominator: scales to 1 amplitude at x=xTap(evalIx)
term = [1 -tapLocations(loopIx+1)] / (tapLocations(evalIx+1)-tapLocations(loopIx+1));
% multiply polynomials => convolve coefficients
polyCoeff = conv(polyCoeff, term);
end
end
% TEST:
% The Lagrange polynomial evaluates to 1 at the location of the tap
% corresponding to evalIx
thisTapLocation = tapLocations(evalIx+1);
pEval = polyval(polyCoeff, thisTapLocation);
assert(max(abs(pEval) - 1) < 1e6*eps);
% The Lagrange polynomial evaluates to 0 at all other tap locations
tapLocations(evalIx+1) = [];
pEval = polyval(polyCoeff, tapLocations);
assert(max(abs(pEval)) < 1e6*eps);
end
# 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')
%For a real-valued, single tone sine or cosine wave, estimate the frequency
%using an fft. Resolution will be limited by the sampling rate of the tone
%and the number of points in the fft.
%Procedure:
%1. Generate the test tone with a given Fs and N
%2. Take the fft and keep the first half
%3. Detect the maximum peak and its index
%4. Equate this index to a frequency
%Error Estimate:
% (Fs/N)= Hz/bin
%Things to keep in mind:
%Adding noise to the input will skew the results
%windowing will help with the fft
%% Clear and close everything
clc;
close all;
clear all; %remove this line if you are trying to use breakpoints
%% Just copy into m file, save and run
Fs = 1e6; %1MHz
fi = 1.333e3; %1kHz
t = 0:1/Fs:.5;
y = sin(2*pi*fi*t);
%Plot the time and frequency domain. Be sure to zoom in to see the waveform
%and spectrum
figure;
subplot(2,1,1);
temp = (2/fi)*Fs;
plot(y);
xlabel('time (s)');
subplot(2,1,2);
sX=fft(y);
N=length(sX);
sXM = abs(sX(1:floor(end/2))).^2; %take the magnitude and only keep 0:Fs/2
plot(0:Fs/N:(Fs/2)-1,sXM);
xlabel('Frequency')
ylabel('Magnitude')
[vv, ii]=max(sXM); %find the index of the largest value in the spectrum
freqEst = (ii-1)*Fs/N; %units are Hz
resMin = Fs/N; %units are Hz
display(freqEst);%display to the command window the results
display(resMin);
# 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 [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)
//Caption:Determination of spectrum of a signal
//With maximum normalized frequency f = 0.1
//using Rectangular window and Blackmann window
clear all;
close;
clc;
N = 61;
cfreq = [0.1 0];
[wft,wfm,fr]=wfir('lp',N,cfreq,'re',0);
disp(wft,'Time domain filter coefficients hd(n)=');
disp(wfm,'Frequency domain filter values Hd(w)=');
WFM_dB = 20*log10(wfm);//Frequency response in dB
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));
end
wft_blmn = wft'.*h_balckmann;
disp(wft_blmn,'Blackmann window based Filter output h(n)=')
wfm_blmn = frmag(wft_blmn,length(fr));
WFM_blmn_dB =20*log10(wfm_blmn);
subplot(2,1,1)
plot2d(fr,WFM_dB)
xgrid(1)
xtitle('Power Spectrum with Rectangular window Filtered M = 61','Frequency in cycles per samples f','Energy density in dB')
subplot(2,1,2)
plot2d(fr,WFM_blmn_dB)
xgrid(1)
xtitle('Power Spectrum with Blackmann window Filtered M = 61','Frequency in cycles per samples f','Energy density in dB')
