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
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);
function [spl, freq] = iso226(phon);
%
% Generates an Equal Loudness Contour as described in ISO 226
%
% Usage: [SPL FREQ] = ISO226(PHON);
%
% PHON is the phon value in dB SPL that you want the equal
% loudness curve to represent. (1phon = 1dB @ 1kHz)
% SPL is the Sound Pressure Level amplitude returned for
% each of the 29 frequencies evaluated by ISO226.
% FREQ is the returned vector of frequencies that ISO226
% evaluates to generate the contour.
%
% Desc: This function will return the equal loudness contour for
% your desired phon level. The frequencies evaulated in this
% function only span from 20Hz - 12.5kHz, and only 29 selective
% frequencies are covered. This is the limitation of the ISO
% standard.
%
% In addition the valid phon range should be 0 - 90 dB SPL.
% Values outside this range do not have experimental values
% and their contours should be treated as inaccurate.
%
% If more samples are required you should be able to easily
% interpolate these values using spline().
%
% Author: sparafucile17 03/01/05
% /---------------------------------------\
%%%%%%%%%%%%%%%%% TABLES FROM ISO226 %%%%%%%%%%%%%%%%%
% \---------------------------------------/
f = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 ...
1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500];
af = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 0.330 0.315 ...
0.301 0.288 0.276 0.267 0.259 0.253 0.250 0.246 0.244 0.243 0.243 ...
0.243 0.242 0.242 0.245 0.254 0.271 0.301];
Lu = [-31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3 -8.1 -6.2 -4.5 -3.1 ...
-2.0 -1.1 -0.4 0.0 0.3 0.5 0.0 -2.7 -4.1 -1.0 1.7 ...
2.5 1.2 -2.1 -7.1 -11.2 -10.7 -3.1];
Tf = [ 78.5 68.7 59.5 51.1 44.0 37.5 31.5 26.5 22.1 17.9 14.4 ...
11.4 8.6 6.2 4.4 3.0 2.2 2.4 3.5 1.7 -1.3 -4.2 ...
-6.0 -5.4 -1.5 6.0 12.6 13.9 12.3];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Error Trapping
if((phon < 0) | (phon > 90))
disp('Phon value out of bounds!')
spl = 0;
freq = 0;
else
%Setup user-defined values for equation
Ln = phon;
%Deriving sound pressure level from loudness level (iso226 sect 4.1)
Af=4.47E-3 * (10.^(0.025*Ln) - 1.15) + (0.4*10.^(((Tf+Lu)/10)-9 )).^af;
Lp=((10./af).*log10(Af)) - Lu + 94;
%Return user data
spl = Lp;
freq = f;
end
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);
