function [b,a] = rc_filter(R, C, Fs, filter_type)
% Returns equivalent IIR coefficients for an analog RC filter
%
% Usage: [B,A] = RC_FILTER(r, c, fs, type);
%
% R is the resistance value (in ohms)
% C is the capacitance value (in farrads)
% FS is the digital sample rate (in Hz)
% type is a character string defining filter type
% Choices are: 'high' or 'low'
%
% Highpass filters have the analog structure:
%
%
% | |
% Vi o--------| |----------+---------o Vo
% | | |
% C |
% ---
% | | R
% | |
% ---
% |
% |
% o---------------------+---------o
% GND
%
%
% Lowpass filters have the analog structure:
%
%
% |-----|
% Vi o--------| |------+---------o Vo
% |-----| |
% R |
% ----- C
% -----
% |
% |
% o---------------------+---------o
% GND
%
% This function uses a pre-calculated equation for both of these circuits
% that only requires the resistance and capacitance value to get a true
% digital filter equivalent to a basic analog filter.
%
% The math behind these equations is based off the basic bilinear transform
% technique that can be found in many DSP textbooks. The reference paper
% for this function was "Conversion of Analog to Digital Transfer
% Functions" by C. Sidney Burrus, page 6.
%
% This is also the equivalent of a 1st order butterworth with a cuttoff
% frequency of Fc = 1/(2*pi*R*C);
%
% Author: sparafucile17 07/01/02
%
% Verify that cutoff of the analog filter is below Nyquist
if( (1/(2*pi*R*C)) > (Fs/2) )
error('This analog filter cannot be realized with this sample rate');
end
% Default to allpass if invalid type is selected
b = [ 1 0 ];
a = [ 1 0 ];
% Constants
RC = R * C;
T = 1 / Fs;
% Analog Cutoff Fc
w = 1 / (RC);
% Prewarped coefficient for Bilinear transform
A = 1 / (tan((w*T) / 2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
% s
% T(s) = -------
% s + 1
%
%
% using Bilinear transform of
%
% 1 ( 1 - z^-1 )
% s --> ----------- * ------------
% tan(w*T/2) ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(strcmp(filter_type,'high'))
b(1) = (A) / (1 + A);
b(2) = -b(1);
a(2) = (1 - A) / (1 + A);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following equations were derived from
%
% 1
% T(s) = -------
% s + 1
%
%
% using Bilinear transform of
%
% 1 ( 1 - z^-1 )
% s --> ----------- * ------------
% tan(w*T/2) ( 1 + z^-1 )
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(strcmp(filter_type,'low'))
b(1) = (1) / (1 + A);
b(2) = b(1);
a(2) = (1 - A) / (1 + A);
end
% *****************************************************************
% Using the Moore-Penrose Pseudoinverse to remove known disturbances from
% measured data
%
% all vectors are column vectors.
% *****************************************************************
close all; clear all;
% *****************************************************************
% Example: A measured signal has been polluted with power-line hum,
% at a fundamental frequency, 3rd and 5th harmonic, with unknown
% amplitude and phase.
% *****************************************************************
m=1000; % number of measured points
fs=22000; % sampling rate
fHum=50; % known 50 Hz power line hum frequency
index=(1:m)';
% create power-line hum terms
hum=0.87654*sin(2*pi*index*fHum/fs+0.12345);
hum3rd=0.65432*sin(2*pi*index*3*fHum/fs+0.45678);
hum5th=0.3321*sin(2*pi*index*5*fHum/fs+0.78901);
% create actual data
idealData=0.1*randn(m, 1);
% create measured data - deteriorated by hum
measuredData=idealData+hum+hum3rd+hum5th;
% *****************************************************************
% Processing to remove the hum
% *****************************************************************
% We know frequency and number of harmonics
% We DO NOT the amplitude and phase of each term.
% note: a*sin(ometaT)+b*cos(omegaT) is equivalent to c*sin(omegaT+d)
% where either [a, b] or [c, d] are the unknown parameters.
basis=[sin(2*pi*index*1*fHum/fs), cos(2*pi*index*1*fHum/fs), ...
sin(2*pi*index*3*fHum/fs), cos(2*pi*index*3*fHum/fs), ...
sin(2*pi*index*5*fHum/fs), cos(2*pi*index*5*fHum/fs)];
% *****************************************************************
% Moore-Penrose Pseudoinverse: Least-squares fit between the basis
% waveforms and the measured signal.
% *****************************************************************
leastSquaresCoeff=pinv(basis)*measuredData;
reconstructedHum=basis*leastSquaresCoeff;
correctedData=measuredData-reconstructedHum;
% *****************************************************************
% Plot measured signal and reconstructed hum
% *****************************************************************
figure(); hold on;
plot(measuredData);
plot(reconstructedHum, 'r', 'LineWidth', 2);
legend({'measured data', 'reconstructed hum signal'});
% *****************************************************************
% Plot the remainder after subtracting the reconstructed hum.
% Compare with ideal result
% *****************************************************************
figure(); hold on;
plot(idealData, 'k', 'LineWidth', 2);
plot(correctedData);
plot(idealData-correctedData, 'r', 'LineWidth', 3);
legend({'ideal result', 'measured and corrected', 'error'})
% David Valencia
% Posted @ DSPRelated.com
% in: http://www.dsprelated.com/showcode/10.php
function [y] = upsample2(x,N)
lx = length(x);
y = 0;
% Put first sample in the first position of the output
y(1) = x(1);
for i=2:lx
% Create a "chunk" (or piece) of vector, composed
% with a series of zeros and the next sample from
% the input
chunk = [zeros(1,N-1),x(i)];
% Append the chunk to the current output vector
y = [y,chunk];
end
function DemoProposed
% Author: J. Selva.
% First version: Sept. 2007.
%
% This m file contains the code for generating Fig. 2 in Sec. IV.A of paper
% "Convolution-based trigonometric interpolation of band-limited signals", IEEE
% Transactions on Signal Processing, vol. 56, n. 11, pp. 5465-5477, Nov. 2008. The
% equation numbers and references in this m file refer to this publication.
%
% Some comments:
%
% -The function 'KnabAPPeriodicApprox' computes the coefficients G(k/NT) in Sec. IV. This
% function takes several minutes, given that it involves a large number of numerical
% integrals which are performed in function KnabAPSpectrum. However, this must be done
% only once for fixed T, B, P and N, given that the G(k/NT) are independent of the sample
% values.
%
% 'RealTestSignal' implements a band-limited signal which consists of a sum of sinc pulses
% with arbitrary amplitudes and delays.
%
% 'FraserFFTUpsampling' implements the well-known zero-padding FFT
% algorithm. 'ProposedFFTUpsampling' implements the upsampling method in Sec. IV. It is
% worth comparing the code of both functions below.
warning off
T = 1; %Sampling period.
B = 0.5; %Two-sided bandwidth.
P = 18; %Pulse semi-length.
N = 1024; %Number of samples.
%Test signal composed of N+200 sinc pulses with random amplitudes and delays.
Sig = @(t)RealTestSignal(t,B,T,[-100,N+100],N+200);
%Compute coefficients Gn. This takes a while (5 minutes) because it is necessary to
%compute a large number of numerical integrals. However, since this step does not depend
%on the signal samples, it can be performed off-line and once ever.
disp(['Computing spectral coefficients. This takes a few minutes (7 or less), but ', ...
'must be done only',sprintf('\n'),'once ever for fixed T, B, P, and N.'])
% 'KnabAPPeriodicApprox' returns the frequency grid, specified by the first frequency fo
% and the frequency increment Df, the coefficients Gn, one for each grid frequency, and
% the truncation index kg.
[fo,Df,Gn,kg] = KnabAPPeriodicApprox(T,B,P,N);
sRef = Sig((0:N-1)*T); %Vector of signal samples with spacing T.
a = 10; %Oversampling factor.
sUpExact = Sig((0:a*N-1)*T/a); %Exact signal samples with oversampling a.
sUpFraser = FraserFFTUpsampling(sRef,a); %Interpolated samples using a zero-padded FFT.
sUpProposed = ProposedFFTUpsampling(sRef,T,B,12,a,Gn,kg); %Interpolated samples using
%the proposed method.
to = 0; %First instant in time grid.
Dt = T/a; %Grid time increment.
%Plot the interpolation error in dB for the zero-padding FFT method and for the proposed
%method.
plot(to+(0:a*N-1)*Dt,20*log10(abs([sUpExact - sUpFraser;sUpExact - sUpProposed])))
xlabel('Time t/T')
ylabel('Interpolation error (dB)')
annotation(gcf,'textarrow',[0.2911 0.2536],[0.8119 0.7429],...
'String',{'This is the zero-padding','FFT error'});
annotation(gcf,'arrow',[0.3321 0.06964],[0.4833 0.6333]);
annotation(gcf,'textarrow',[0.3446 0.2643],[0.419 0.3548],...
'String',{'This is the error of the','proposed FFT method. Note',...
'that the scale is in dB !!'});
warning on
%=====================================================================
%
% Conventional zero-padding FFT algorithm
function sUp = FraserFFTUpsampling(s,a)
N = length(s);
sF = a * fftshift(fft(s)); %Compute the FFT of the signal samples.
if mod(N,2) == 0 %Pad with zeros to obtain over-sampling factor a.
sF(end) = 0.5*sF(end);
sF = [zeros(1,N*(a-1)/2),sF,zeros(1,N*(a-1)/2)];
else
sF = [zeros(1,(N*(a-1)+1)/2),sF,zeros(1,(N*(a-1)-1)/2)];
end
sUp = ifft(ifftshift(sF)); %Compute the inverse FFT.
%=====================================================================
%
% Proposed FFT algorithm
function sUp = ProposedFFTUpsampling(s,T,B,P,a,Gn,kg)
N = length(s);
sF = fft(s); %Compute the FFT of the signal samples.
sF = sF(1+mod(-kg:kg,N)) .* Gn * (N*a); %Repeat some of the spectral samples periodically,
%and apply the window specified by the vector Gn.
sF = sF(:).';
if mod(N,2) == 0 %Pad with zeros to obtain over-sampling factor a.
sF = [zeros(1,N*a/2-kg),sF,zeros(1,N*a/2-kg-1)];
else
sF = [zeros(1,(N*a-2*kg-1)/2),sF,zeros(1,(N*a-2*kg-1)/2)];
end
sUp = ifft(ifftshift(sF)); %Compute the inverse FFT.
%=====================================================================
%
% Returns the spectral coefficients Gn and the frequency grid specification f0, Df, Gn.
function [f0,Df,Gn,n0] = KnabAPPeriodicApprox(T,B,P,N)
Df = 1/(N*T);
n0 = floor((1/T-B/2)/Df);
f0 = -Df*n0;
Bg = 2/T-B;
Gn = KnabAPSpectrum(f0+Df*(0:2*n0),B,Bg,(P+1/2)*T)/N;
%=====================================================================
%
% Auxiliary function that computes the spectral coefficients Gn.
function G = KnabAPSpectrum(f,B,Bg,TL)
% B: Two-sided signal bandwidth.
% Bg: Two-sided bandwidth of output pulse.
% Be: Two-sided bandwidth of fast-decay pulse.
% TL: Window's width between zeros in time domain.
%
% The two-sided bandwidth of the output pulse is Bg = B+2*Be.
Nf = numel(f);
G = zeros(size(f));
Be = (Bg-B)/2;
Ba = (Bg+B)/2;
Int = @(x)ConvWindow(x,Be,2*TL);
for k = 1:Nf
if f(k)-Ba/2 <= -Be/2 && Be/2 <= f(k)+Ba/2
G(k) = 1;
else
IInt = IntervalIntersection([-Be/2,Be/2],f(k)+[-Ba/2,Ba/2]);
if ~isempty(IInt)
G(k) = quad(Int,IInt(1),IInt(2),1e-16);
end
end
end
% Auxiliary function. It computes samples of the continuous convolution of a square pulse
% with the Kaiser window.
function W = ConvWindow(f,B,TL)
%Kaiser window.
%B: Two-sided window width.
%TL: Width between zeros in time domain.
W = zeros(size(f));
I = abs(f) < B/2;
W(I) = besseli(0,pi*B*(TL/2)*sqrt(1-(2*f(I)/B).^2))/(B*sinc(j*B*TL/2));
function [IInt,RelAInd,RelBInd] = IntervalIntersection(IA,IB)
if IA(2) < IB(1) || IB(2) < IA(1) %No overlap case.
IInt = [];
RelAInd = [];
RelBInd = [];
return
end
SwapIntervals = diff(IA) > diff(IB);
if SwapIntervals %Force IA to be the sorter interval.
Ipr = IA;
IA = IB;
IB = Ipr;
end
if IA(1) < IB(1)
IInt = [IB(1),IA(2)]; %There is intersection to the left.
RelAInd = IB(1)-IA(1);
RelBInd = 0;
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
return
end
if IA(2) <= IB(2)
IInt = IA; %IA lies inside IB.
RelAInd = 0;
RelBInd = IA(1)-IB(1);
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
return
end
IInt = [IA(1),IB(2)]; %There is intersection to the right.
RelAInd = 0;
RelBInd = IA(1)-IB(1);
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
%=====================================================================
%
% Computes samples of a signal formed by sincs with random amplitudes and delays.
function s = RealTestSignal(t,B,T,I,Np)
Sr = rand('seed');
rand('seed',0); %Always select the same random numbers in the sequel.
M = numel(t);
Ap = (rand(Np,1)-0.5)*2 ; %Random amplitudes.
tp = I(1)+rand(Np,1)*diff(I); %Random delays.
rand('seed',Sr)
st = size(t);
t = t(:);
s = zeros(numel(t),1);
for r = 1:Np
s = s + sinc(B*(t-tp(r)))*Ap(r); %Compute sum of sincs
end
s = reshape(s,st);
function [S] = dct_2d(P) %where P is the matrix to apply the dct
N=length(P);
Cu=[1/sqrt(2),ones(1,N-1)]; %Constant coeficients d1
Cv=[1/sqrt(2),ones(1,N-1)]; %Constant coeficients d2
S=zeros(N,N);
for v=0: N-1;
for y=0: N-1;
for u=0: N-1;
for x=0: N-1; %serial processing
S(u+1,v+1)=S(u+1,v+1)+Cu(u+1)*Cv(v+1)*P(x+1,y+1)*cos((2*x+1)*pi*u/(2*N))*cos((2*y+1)*pi*v/(2*N));
end;
end;
end;
end;
S=(2/N)*S; %final result
% Usage: [B,A] = COMB(order, scalar);
%
% ORDER is the number of samples delayed prior to add
% SCALAR is the coefficient that will be applied to
% the delayed signal path at the final summation block.
%
% Note, there are two types of comb filters. A DC-blocker and a DC-passer.
% To get a DC-Blocker (tooth at DC), pass in a -1 for the scalar.
% To get a DC-Passer (+6dB at DC), pass in a +1 for the scalar.
%
% By default, if the scalar is not passed, a DC-Passer is assumed.
%
% Author: sparafucile17 03/16/04
% Validate that the proper argument count was supplied
error(nargchk(1, 2, nargin));
% Use scalar is passed as an argument, otherwise assume scalar=1;
if (nargin == 1)
scalar = 1;
else
scalar = varargin{1};
end
% Input has zeros to simulate a single delay of N samples
a = [ 1 ];
b = [ 1 zeros(1, N-1) scalar*1];
function Demo
%Author: J. Selva.
%Date: June 2011.
close all
format short g
format compact
T = 1; %Sampling period in both dimensions.
B = 1/(1.223*T); %This is the two-sided bandwidth in both dimensions. It is a typical ...
%value in SAR image co-registration.
P = 18; %Kernel semi-length.
Q = 10; %Polynomial order in Farrow structure.
%The next call to "ChebApproxPolys" computes all Farrow coefficients.
t0 = -P*T-T/2; %The parameters t0, Nt, 2*P+1 specify the 2*P+1 intervals in which Knab's ...
%pulse is to be approximated by polynomials.
Nt = 2*P+1;
Dt = T;
IOut = [-T/2,T/2];
%Perform Chebyshev interpolation. The parameters T, B and P are passed to KnabAPPulse.
K = ChebApproxPolys('KnabAPPulse',t0,Nt,Dt,IOut,40,Q,T,B,P).';
%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = 100;
x0 = (0:Lx-1)';
x0 = x0(:,ones(1,Lx));
y0 = x0.';
x0 = x0(:);
y0 = y0(:);
disp('Computing image in a regular grid...')
tic;
Image = reshape(TestImage5(x0,y0,B),[Lx,Ly]);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
clear x0 y0
%Specify NP points on the image at random.
NP = 10000;
x = P*T+rand(1,NP)*(Lx-2*P)*T;
y = P*T+rand(1,NP)*(Ly-2*P)*T;
%Evaluate image directly.
disp('Computing image at selected positions exactly...')
tic;
vRef = TestImage5(x,y,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
%Interpolate the image using the 2D-Farrow structure.
disp('Computing image at selected positions using Farrow 2-D...')
tic;
v = Farrow2DInterp(Image,K,x,y);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
%Show the error.
disp(['Maximum interpolation error: ', ...
num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])
disp(' ')
disp('Display image and interpolation error:')
disp(' ')
%Evaluate the test image in a regular grid with spacing 1.
Lx = 100;
Ly = Lx;
x1 = (P*T:0.5*T:(Lx-P)*T)';
x1 = x1(:,ones(1,length(x1)));
y1 = x1.';
x1 = x1(:);
y1 = y1(:);
tic;
disp('Computing image in over-sampled grid ...')
vRef = TestImage5(x1,y1,B);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
tic;
disp('Computing signal in over-sampled grid using 2-D Farrow ...')
v = Farrow2DInterp(Image,K,x1,y1);
disp(['(Timing: ',num2str(toc), ' sec.)'])
disp(' ')
x1 = (P*T:0.5*T:(Lx-P)*T)';
y1 = x1.';
disp(['Maximum interpolation error: ', ...
num2str(20*log10(max(abs(vRef-v))/max(abs(vRef)))),' dB'])
vRef = reshape(vRef,[length(x1),length(y1)]);
v = reshape(v,[length(x1),length(y1)]);
mesh(x1,y1,abs(vRef))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Image''s abs. value')
title('Test Image (abs. value)')
figure
mesh(x1,y1,20*log10(abs(v-vRef)/max(abs(vRef(:)))))
xlabel('x coordinate')
ylabel('y coordinate')
zlabel('Interp. error (dB)')
title('Interpolation error (dB)')
function P = ChebApproxPolys(Func,t0,Nt,Dt,IOut,NPoints,Q,varargin)
%This function constructs Chebyshev interpolation polynomials in a set of consecutive ...
%intervals.
%
% Func is the function to be interpolated.
% t0 is the first instant of the first interval.
% Nt is the number of intervals.
% IOut is the interval in which the output polynomials will be evaluated.
% Dt is the interval width. The intervals are [t0,t0+Dt], [t0+Dt,t0+2*Dt],....
% [t0+(Nt-1)*Dt,t0+Nt*Dt].
% NPoints is the number of points in which Func is evaluated, for each interval.
% Q is the output polynomial order.
% varargin contains extra parameters which are passed to varargin.
%Construct matrix that maps function values to polynomial coefficients.
%The algorithm is that in section 9.8 of
%
% W. H. Press et al., Numerical Recipes in C. Cambridge University
% Press, 1997.
M = MapPolyMatrix([-1,1],IOut,Q) * ChebCoefMatrix(Q) * ...
BlockChebInvMatrix(NPoints,Q);
% Get Chebyshev roots.
r = ChebRoots(NPoints);
r = r(:);
P = zeros(Q,Nt);
%Compute polynomial for each interval.
for k = 0:Nt-1,
v1 = feval(Func,MapIaToIb(r,[-1,1],t0+[k,k+1]*Dt),varargin{:});
P(:,k+1) = double(M * v1(:));
end
function M = MapPolyMatrix(IIn,IOut,Q)
%For a polynomial with Q coefficients which is evaluated in interval IIn, this function ...
%gives the matrix that transforms its coefficients, so that it can be evaluated in the ...
%associated points in interval IOut.
ayx = (IIn(1)-IIn(2))/(IOut(1)-IOut(2));
byx = IIn(1)-ayx*IOut(1);
if byx == 0
M = diag(ayx.^(0:Q-1));
else
M = zeros(Q);
for q = 0:Q-1
for r = 0:q
M(r+1,q+1) = nchoosek(q,r) * ayx^r * byx^(q-r);
end
end
end
function M = ChebCoefMatrix(N)
M = zeros(N);
M(1) = 1;
if N > 1
M(2,2) = 1;
end
for n = 3:N,
M(2:n,n) = 2*M(1:n-1,n-1);
M(1:n-2,n) = M(1:n-2,n) - M(1:n-2,n-2);
end
function M = BlockChebInvMatrix(N,Q)
r = ChebRoots(N);
r= r(:);
M = zeros(N,Q);
if Q >=1
M(:,1) = 1;
end
if Q >= 2
M(:,2) = r;
end
for q = 3:Q,
M(:,q) = 2 * r .* M(:,q-1) - M(:,q-2);
end
M = M.'/N;
M(2:end,:) = M(2:end,:) * 2;
function R = ChebRoots(N)
R = cos(pi*((1:N)-1/2)/N);
function y = MapIaToIb(x,Ia,Ib)
%Linear transformation that maps interval Ia onto interval Ib.
y = (Ib(1)+Ib(2))/2 + (x-(Ia(1)+Ia(2))/2) * ((Ib(2)-Ib(1))/(Ia(2)-Ia(1)));
function [g,BoundAP] = KnabAPPulse(t,T,B,P)
%Knab interpolation pulse in
%
% J. J. Knab, 'Interpolation of band-limited functions using the Approximate
% Prolate series', IEEE Transactions on Information Theory, vol. IT-25,
% no. 6, pp. 717-720, Nov 1979.
%
%Grid points -P*T:T:P*T.
%Interpolation interval: [-T/2,T/2].
g = sinc(t/T) .* sinc((1/T-B)*sqrt(t.^2-(P*T)^2)) / sinc(j*(1-B*T)*P);
if nargout == 2
BoundAP = 1/sinh(P*pi*(1-B*T));
end
function v = TestImage5(x,y,B)
%Test image formed by a sum of random sincs.
Seed = rand('state');
rand('state',362436069);
L = 100; %Image size in both dimensions.
Np = 5000;
p = rand(Np,2)*L;
A = exp(j*2*pi*rand(Np,1));
rand('state',Seed);
v = zeros(size(x));
for k = 1:length(x(:))
v(k) = A(:).' * (sinc(B*(x(k)-p(:,1))) .* sinc(B*(y(k)-p(:,2))));
end
function vx = Farrow2DInterp(Image,Kernel,x,y)
%Farrow interpolation in two dimensions using interlaced FFTs.
vxSize = size(x);
[LK,NK] = size(Kernel);
[Lx,Ly] = size(Image);
nx = round(x);
ny = round(y);
ux = x-nx;
uy = y-ny;
np = 1+nx+ny*Lx;
np = np(:);
clear x y nx ny
P = (LK-1)/2;
NFFTx = 2^ceil(log2(Lx+LK-1));
NFFTy = 2^ceil(log2(Ly+LK-1));
KernelFx = fft(Kernel,NFFTx);
KernelFy = fft(Kernel,NFFTy);
ImageFx = fft(Image,NFFTx);
vx = zeros(size(np))+j*zeros(size(np));
for kx = NK:-1:1,
ImageCx = ifft(ImageFx .* KernelFx(:,kx(ones(1,Ly))),NFFTx);
ImageCx = ImageCx(1+P:Lx+LK-1-P,:);
ImageCxFy = fft(ImageCx,NFFTy,2);
vy = zeros(size(np))+j*zeros(size(np));
for ky = NK:-1:1,
ImageCxCy = ifft(ImageCxFy .* KernelFy(:,ky(ones(1,Lx))).',NFFTy,2);
ImageCxCy = ImageCxCy(:,1+P:Ly+LK-1-P);
vy(:) = vy(:) .* uy(:) + ImageCxCy(np);
end
vx(:) = vx(:) .* ux(:) + vy(:);
end
vx = reshape(vx,vxSize);
clf;
% First plot
t = 0:0.0005:1;
f = 5;
x_t = cos(2*pi*f*t);
subplot (4,1,1);
plot (t,x_t);
grid;
xlabel ('Time [sec]');
ylabel ('Amplitude');
title ('Continuous-time signal x(t)');
axis([0 1 -1.2 1.2])
% Second plot
n = 0:0.125:1;
x_n = cos(2*pi*f*n);
subplot(4,1,2);
y = zeros(1,length(t));
for i = 1:length(n)
y = y + x_n(i)*sinc(t/0.125 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs < bound] (T=0.125)');
axis ([0 1 -1.2 1.2]);
% Third plot
n = 0:0.1:1;
x_n = cos(2*pi*f*n);
subplot(4,1,3);
y = zeros(1,length(t));
for i = 1:length(n)
y = y + x_n(i)*sinc(t/0.1 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs = bound] (T=0.1)');
axis ([0 1 -1.2 1.2]);
% Fourth plot
n = 0:0.075:1;
x_n = cos(2*pi*f*n);
subplot(4,1,4);
y = zeros(1,length(t));
for i = 1:length(n)
y = y + x_n(i)*sinc(t/0.075 - i + 1);
end
plot(t,y);
grid;
xlabel ('Time, sec');
ylabel ('Amplitude');
title ('Reconstructed continuous-time signal y(t) [fs > bound] (T=0.075)');
axis ([0 1 -1.2 1.2]);
function [] = SineSweep(fstart, fstop, length, method, fs, bits, PHI)
%SINESWEEP Creates a custom sine wave sweep.
% SineSweep(fstart, fstop, length, [method], [fs], [bits], [PHI])
% fstart (Hz) - instantaneous frequency at time 0
% fstop (Hz) - instantaneous frequency at time length
% length (s) - the length of time to perform the sweep
% method - 'quadratic', 'logarithmic', or 'linear' (default)
% fs (Hz) - the sampling frequency (default = 48KHz)
% bits - bit depth of each sample 8, 16 (default), 24, or 32
% PHI (deg) - initial phase angle. cos = 0, sin = 270 (default)
%
% Written by: Ron Elbaz
% Date: February 15, 2011
%check and set missing parameters
if exist('method','var') ~= 1
method = 'linear';
end
if exist('fs','var') ~= 1
fs = 48000;
end
if exist('bits','var') ~= 1
bits = 16;
end
if bits ~= 8 && bits ~= 16 && bits ~= 24 && bits ~= 32
bits = 16;
end
if exist('PHI','var') ~= 1
PHI = 270;
end
%make sure start and stop frequencies are valid
if fstop < fstart
temp = fstart;
fstart = fstop;
fstop = temp;
end
%avoid aliasing
if fstart > fs/2
fstart = fs/2;
end
if fstop > fs/2
fstop = fs/2;
end
%create time vector
t = 0:(1/fs):length;
%scale to avoid clipping due to rounding error
Y = 0.9999*chirp(t,fstart,length,fstop, method, PHI);
%play the sound
% sound(Y,fs);
%store sweep to wav file
wavwrite(Y,fs,bits,'SineSweep.wav');
end
%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
%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
%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);
