Sampling
The traditional approach to sampling theory starts by considering the circumstances under which a continuous signal can be reconstructed from knowledge of its values only at certain sample points, equally spaced in time (or space, or whatever).
The sampling process is modeling by multiplying the original signal by an impulse train, which is another continuous signal that has the value zero except at a set of discrete, equally-spaced points. We then consider the (continuous) Fourier transform of the result, and this leads directly to the sampling theorem, which provides a simple and clear definition of the circumstances under which the original signal can be recovered, and also give a simple prescription for the recovery.
In the digital domain, someone has usually taken care of "sampling" -- in this sense -- already. Old-fashioned sampling theory might come up if we need to decide what regular sampling frequency we should use to deal with some real-world analog signal, or how to deal with an analog signal that's been sampled in an unusual way, say at irregular intervals. But most of today's samping problems involve digital transformations of digital data. One common case is sample rate changing -- how do we change a regularly-sampled signal from some sampling frequency to another?
In these notes, we'll consider the problems of sampling in two other ways, both entirely in the discrete domain. One discussion explores the relationship between sampling issues and modular arithmetic, as seen in the case of the wavetable oscillator; and the other discussion explores the relationship between sampling and projection into and out of a subspace.
Resampling and the wavetable oscillator.
Suppose that we have a vector P of length N, containing whatever.
If we repeat P over and over again, forming an indefinitely long vector called (say) Pstar, it's obvious that the result is periodic, with period N. That is , if k = j + n*N, for any integer n, then Pstar(k) == Pstar(j).
P might contain one period of a sine wave, or one period of a violin sound, or one period of a gorilla grunt, or just N random numbers---Pstar is still periodic with period N.
Since each period of Pstar is exactly the same, by construction, Pstar(j) = P(modulo(j,N)).
N.B. To make this true in MATLAB we would have to use a special modulo() function whose domain is from 1 to N rather than from 0 to N-1.
If we regard the vector P as a digital version of a sound, what is the frequency of the periodic version Pstar? To answer this question, we have to know the time interval Tsamp between the digital samples, or equivalently the sampling frequency Sfreq = 1/Tsamp.
If Tsamp is 1/10000 seconds (Sfreq = 10000 Hertz) then one period of Pstar is N/10000 seconds long, and Pstar has a frequency of 10000/N.
In the discussion that follows, we will mostly be thinking of the the sampling frequency as a fact about the process of "ditigal to analog conversion" (D/A), i.e. the process of playing the sound. It might equally apply to the process of analog-to-digital conversion (A/D), and the two rates need not be the same.
Thus we could record a violin at 100,000 samples per second, and in this case a note played at 400 Hertz will have a period of 100000/400 = 250 samples, so that our vector P will have length N = 250. If we play Pstar back at 10,000 samples per second, then its frequency will be 10000/250 = 40 Hertz. Here no digital processing needs to be done to change the sampling rate (and hence the frequency content of the signal), we are just running the A/D process by different rules from the D/A process, producing the sound of an enormous violin.
If P is entirely synthetic, then its "natural" period is purely a matter of convention, and the period of the sound corresponding to Pstar will just depend on the output sampling rate. Suppose P contains N samples of a continuous function that we might have chosen to sample at any arbitrary density, so that our vector representing one period might have contained any other number of samples M.
This function might be something inherently periodic like sin(), or it might be something like y = -2x + 1, 0 <= x < 1. In this last example, Pstar will be a sawtooth:
N = 20; x1 = (0:(N-1))/N; y1 = -2*x1 + 1; y1star = [y1 y1 y1 y1]; plot(1:(N*4),y1star,'g-', 1:(N*4),y1star,'bo'); M = 10; x2 = (0:(M-1))/M; y2 = -2*x2 + 1; y2star = [y2 y2 y2 y2]; figure(2); newplot; plot(1:(M*4),y2star,'g-', 1:(M*4),y2star,'bo');
If the playback sampling rate is 10 KHz, then these two examples will have frequencies 10000/20 = 500 Hz and 10000/10 = 1000 Hz.
Because of the way that we set the example up, y2star is identical a version of y1star with every other sample left out:
>> y1star(1:2:80) - y2star
ans =
Columns 1 through 12
0 0 0 0 0 0 0 0 0 0 0 0
Columns 13 through 24
0 0 0 0 0 0 0 0 0 0 0 0
Columns 25 through 36
0 0 0 0 0 0 0 0 0 0 0 0
Columns 37 through 40
0 0 0 0
In general, if we have a table P of N equally-spaced samples of some function, then we can create a table of M equally-spaced samples, for M = N/k, by pulling out every kth sample.
Thus from a table of 1200 samples, we can create tables of 600, 400, 300, 240, 200 etc. for k = 2, 3, 4, 5, 6
If our playback sampling frequency is S, these correspond to pitches of S/1200, S/600, S/400, S/300, S/240 etc.
This process produces exactly what we would have gotten by creating these tables over again from first principles each time, and will be computationally convenient if the function is somewhat slow to compute, like sin() -- the wavetable lookup can be implemented efficiently in a loop that successively increments an index by k, checking against N to see when "wrap-around" is needed.
In cases where the "function" is not known, but only empirically observed -- e.g. the sound signal corresponding to a piano, or a violin, or an oboe -- this method of producing new sampled signals is especially attractive. A variant of this process has been used in many commercial synthesizers, sound cards and synthesis software.
There is one big problem -- the frequencies available from integer values of k are only a few of the ones that we want.
However, there is no reason in principle not to make k -- which we can think of as the skip factor in Pstar -- something other than an integer.
Here's a MATLAB function that calculates k for (almost-) arbitrary resampling. It is set up so that we can think of it as an oscillator driven by a table containing one period of the desired function, and producing n samples at any frequency f for any output sampling rate r:
function s = wosc(P,f,r,n)
%
% Wavetable Oscillator
%
% s = wosc(P,f,r,n)
%
% P = one period of a sampled signal
% f = desired output frequency
% r = output sampling rate
% n = number of output samples to compute
%
% This version just takes the closest value with no interpolation.
lP = length(P);
inc = lP*f/r % increment in P per output sample:
% each output sample is 1/r long;
% each sample in P is nominally (1/f)/lP long;
% thus if intersample distance in P is 1,
% intersample distance in output is (1/r)/((1/f)/lP)
ind = mod(round(0:inc:(inc*n-1)), lP) + 1;
s = P(ind);
return
Let's look at how it works in some cases where k is an integer. Suppose we use a table of 1200 samples of one period of a sine wave, and then use wosc to generate 1200 samples of a sine wave with frequency 10 at sampling rate 12000. This corresponds to k = 1, and should obviously work perfectly if we compare it to what we would have gotten by using the sin() function instead:
P = makesin(1200,1,1,0);
wP = wosc(P,10,12000,1200);
sqrt(sum((P-wP).^2)/1200)
ans =
0
Now let's try to use the same table to generate 1200 samples of a sine wave with frequency 20 at sampling rate 12000. This makes k = 2 and should also work perfectly:
P2 = makesin(1200,2,1,0); wP2 = wosc(P,20,12000,1200); sqrt(sum((P2-wP2).^2)/1200) ans = 3.6622e-16
The reason that it's not completely perfect is that the Matlab sin() function is having to calculate an output for inputs above 2*pi -- probably by projecting them back into the range from 0 to 2*pi -- and this opens up the opportunity for error to creep in due to limited-precision arithmetic. Still, the difference is very small -- and in this case, the wavetable oscillator is probably more accurate!
Now let's try something more significant: we'll calculate the frequency of middle C, and create 12,000 samples of a middle-C sine wave at sampling frequency 12,000 (one second of a middle C sine tone:
>> semitone = 2^(1/12);
>> A = 440;
>> C = (A/2)*semitone^3
C =
261.6256
>> P3 = makesin(12000,C,1,0);
>> wP3 = wosc(P,C,12000,12000);
inc =
26.1626
>> newplot; plot(P3(1:50)); hold on; plot(wP3(1:50),'--'); plot(P3(1:50)-wP3(1:50),'o');
>> title('First 50 samples comparing wosc() with sin()'); hold off;
>> newplot; hold on; plot(P3(11000:11050)); plot(wP3(11000:11050),'--');
>> plot(P3(11000:11050)-wP3(11000:11050),'o');
>> title('Comparing 50 samples starting at index 11,000'); hold off;
In the two plots below, the circles plotted near the zero line represent the differences between the two versions of "middle C." Graphically, the difference does not look like much:
However, there is really a difference between the wosc() and sin() versions of middle C:
>> sqrt(sum((P3-wP3).^2)/12000)
ans =
0.0011
>> max(abs(P3-wP3))
ans =
0.0026
How important is this? Well, .0026 at the scale of +-1 would be about .7 at the scale of +-128, as given by the 8-bit audio that used to be standard in personal computer audio systems. This means that at most the low-order bit will be wrong, which probably doesn't matter very much. On the other hand, .0026 is about 86 at the scale of +-2^15, which is the range used by CD audio and most contemporary digital sound systems. This means that the low-order 6 or 7 bits might be messed up, and this probably would matter, depending on the acoustic material, the listening conditions, and the quality of the audio system being used.
We could get less distortion by using a bigger wavetable. Alternatively, we could do some kind of interpolation, so that rather than just taking the closest (in time) available sample value, we invented a new sample value based on the adjacent samples:
Here is a version of wosc() that does linear interpolation:
function s = iwosc(P,f,r,n)
%
% Interpolating Wavetable Oscillator
%
% s = wosc(P,f,r,n)
%
% P = one period of a sampled signal
% f = desired output frequency
% r = output sampling rate
% n = number of output samples to compute
%
% This version does linear interpolation
lP = length(P);
P1 = [P P(1)]; % in case of interpolation across the wrap
inc = lP*f/r % increment in P per output sample:
% each output sample is 1/r long;
% each sample in P is nominally (1/f)/lP long;
% thus if intersample distance in P is 1,
% intersample distance in output is (1/r)/((1/f)/lP)
ind = mod(0:inc:(inc*n-1), lP) + 1;
ind1 = floor(ind);
ind2 = ceil(ind);
frac2 = ind - ind1;
frac1 = 1 - frac2;
s = frac1.*P1(ind1) + frac2.*P1(ind2);
return
Now let's try the same thing over again:
>> wiP3 = iwosc(P,C,12000,12000); inc = 26.1626 >> sqrt(sum((P3-wiP3).^2)/12000) ans = 1.7694e-06 >> max(abs(P3-wiP3)) ans = 3.4266e-06
This is giving us a maximum error of around 3.4 parts in a million. In terms of 16-bit audio, the maximum error would be about .23, i.e. there will almost never be any difference, and at most it will be in one low order bit. This will be entirely inaudible under nearly all circumstances.
Nevertheless, this type of interpolation is Not Right.
Before we go on to see why, and to see what the Right Answer is, let's consider a more important concern about wavetable oscillators.
Let's go back to the case of a table containing N samples of one period of a sine wave. For integer values of the skip factor k, we decided that the wavetable oscillator was doing exactly the right thing. It should be obvious that there are some limits to this generalization. For instance, if the skip factor is exactly N, then there will be no oscillation at all -- the same sample will be chosen in every period! This is probably not what we wanted.
Since the formula for the increment is k = N*f/S where N is table length f is desired output frequency S is output sampling rate the skip factor will be exactly N where f == S.
We will have exactly the same problem where f == n*S, for any integer n, since this will result in a skip factor of n*N, and for any integer value of n this is the same as a skip factor of 0. Notice that this has nothing to do with the size of the wavetable -- making the wavetable longer or shorter doesn't help or hurt.
We can also get into trouble for other skip factors.
Consider the case where two skip factors differ by n*N/2. Thus we might have N = 40, S = 100, and f1 = 25, f2 = 75. This produces k1 = 40*25/100 = 10, or k2 = 40*75/100 = 30. The start of the index sequence in the two cases will be: x = 0:10; s1 = mod(10*x,40) indices for frequency 25 (+ 1 for MATLAB!) s2 = mod(30*x,40) indices for frequency 75 (+ 1 for MATLAB!)
s1 looks OK, since its period is five, and at sampling frequency 100, a period of 5 corresponds to a frequency of 25, which is what we asked for.
However, something is wrong with s2. Its period is also 5! so it repeats 25 times in 100 samples, not 75 times. Nevertheless, s2 is not exactly the same as s1 -- the indices run in reverse order, so that the resulting oscillation is a time-reversed version of the oscillation at frequency 25. Therefore the result will correspond to an oscillation of our wavetable's function at a frequency of -25.
This same kind of thing -- aliasing with a (possibly time-reversed) version of a different frequency -- will happen any time we ask for a frequency that translates to a skip factor greater than N/2.
For instance, suppose we have N = 40, S = 20, f1 = 2, f2 = 18. This produces k1 = 40*2/20 -> 4, k2 = 40*18/20 -> 36 x = 0:10; s1 = mod(4*x,10) indices for frequency 2 s2 = mod(36*x,10) indices for frequency 18
This time around, frequency 18 has been aliased into frequency -2.
This is a simple consequence of the nature of modular arithmetic, and modular arithmetic in turn comes into play here because we are sampling a PERIODIC signal, so that indices are equivalent modulo the signal's period. The notion of periodicity in turn is inevitable once we start talking about FREQUENCIES.
We observed that aliasing will occur whenever we get a skip factor greater than N/2, i.e. N*f/S > N/2 where f is the desired output frequency S is the sampling rate N is the number of samples in the wave table A bit of algebra turns this into f > S/2
Thus, if we try to use a wavetable oscillator to create a frequency f at a sampling rate S, what we actually get is a frequency
fhat = (f mod S) if (f mod S) < S/2
S/2 - (f mod S) if (f mod S) > S/2
For (f mod S) == S/2, either definition applies.
This form of the generalization has nothing special to do with the wavetable oscillator algorithm -- it is a general fact about the representation of periodicity in any sampled signal, no matter how it come about. We've seen it before, in the derivation of the negative frequencies in the DFT. Therefore it is worth memorizing the general case.
Interpolation and decimation in Matlab
Octave (and Matlab) provide us with functions that are intended to change sampling rates of signals in the Right Way, called interp() (for "interpolate") and decimate(). Let's experiment briefly with interp() and decimate().
help interp
INTERP Resample data at a higher rate using lowpass interpolation.
Y = INTERP(X,R) resamples the sequence in vector X at R times
the original sample rate. The resulting resampled vector Y is
R times longer, LENGTH(Y) = R*LENGTH(X).
A symmetric filter, B, allows the original data to pass through
unchanged and interpolates between so that the mean square error
between them and their ideal values is minimized.
Y = INTERP(X,R,L,ALPHA) allows specification of arguments
L and ALPHA which otherwise default to 4 and .5 respectively.
2*L is the number of original sample values used to perform the
interpolation. Ideally L should be less than or equal to 10.
The length of B is 2*L*R+1. The signal is assumed to be band
limited with cutoff frequency 0 < ALPHA <= 1.0.
[Y,B] = INTERP(X,R,L,ALPHA) returns the coefficients of the
interpolation filter B. See also DECIMATE.
help decimate
DECIMATE Resample data at a lower rate after lowpass filtering.
Y = DECIMATE(X,R) resamples the sequence in vector X at 1/R times the
original sample rate. The resulting resampled vector Y is R times shorter,
LENGTH(Y) = LENGTH(X)/R.
DECIMATE filters the data with an eighth order Chebyshev type I lowpass
filter with cutoff frequency .8*(Fs/2)/R, before resampling.
Y = DECIMATE(X,R,N) uses an N'th order Chebyshev filter.
Y = DECIMATE(X,R,'FIR') uses the 30 point FIR filter generated by
FIR1(30,1/R) to filter the data.
Y = DECIMATE(X,R,N,'FIR') uses the N-point FIR filter.
Note: For large R, the Chebyshev filter design might be incorrect
due to numeric precision limitations. In this case, DECIMATE will
use a lower filter order. For better anti-aliasing performance, try
breaking R up into its factors and calling DECIMATE several times.
See also INTERP, RESAMPLE.
Now let's get some signals ready to work with:
s_100_1 = makesin(100,1,1,0); s_100_5 = makesin(100,5,1,0); s_100_19 = makesin(100,19,1,0); s_100_25 = makesin(100,25,1,0); s_20_1 = makesin(20,1,1,0); s_20_5 = makesin(20,5,1,0); s_500_1 = makesin(500,1,1,0); s_500_25 = makesin(500,25,1,0); % ind is the set of index values that pick every 5th sample... ind = 1:5:100;
Just plain desampling a frequency-1 sine wave, by taking every 5th sample, works fine!
>> x = s_100_1(ind);
>> max(abs(x-s_20_1))
ans =
0
MATLAB's decimate function actually screws up here, due to phase shift and edge effects in its fancy Chebyshev filter:
>> x = decimate(s_100_1,5);
>> max(abs(x-s_20_1))
ans =
0.2481
>> newplot; hold on; plot(x,'-'); plot(s_20_1,':')
>> title('Correct sine wave (dots) vs. MATLAB decimate version');
We can do a little better by asking MATLAB to use a Finite Impulse Response filter (which is actually the default for Octave):
x = decimate(s_100_1,5,10,'FIR'); max(abs(x - s_20_1)) newplot; plot(x);
Now let's add up a frequency-1 and a frequency-19 sine wave:
newplot; ss = s_100_1 + s_100_19; plot(ss);
Oops -- if we desample by 5, suddenly there's nothing left!
newplot; plot(ss(ind)); axis([1 20 -1 1]);
Here's why -- the desampled version of frequency 19 is actually aliased to -1 -- and frequency 1 and frequency -1 sine waves cancel:
newplot; plot(s_100_19(ind));
title('s_100_19 desampled by 5 -- turns into sin(-x)');
newplot; plot(s_100_1(ind));
title('s_100_1 desampled by 5 -- still sin(x)');
Now the reason for the low-pass filtering becomes clear: first we need to get rid of frequencies that might alias -- those above S/2 where S is the output sampling rate... then we can desample safely. If we do this, the frequency-19 sine wave is just gone, instead of being aliased.
figure(1);
r2 = decimate(ss,5,10,'FIR');
newplot; plot(r2);hold on; plot(s_20_1,':');
title('Sin_1 plus sin_19, filtered before decimation');
We'll have more to say about upsampling later on. For now, here's a picture of the filter MATLAB uses to interpolate on upsampling.
[y,b] = interp(s_20_1,5); newplot; lplot(b);
A subspace account of sampling
We generate 9 random basis vectors of length 32.
basis = rand(32,9);
Although it is not very interesting, we can plot them:
figure(1); plot(basis);
They look suitably random.
Now we construct a signal of length 32 by taking a random linear combination of our nine basis vectors.
coeffs = rand(9,1) signal = basis*coeffs;
Let's plot the result:
plot(1:32,signal,'-', 1:32,signal,'o');
This signal also looks suitably random.
However, as we know, it is not based on 32 random choices, but only on 9 --- the 9 elements of the vector 'coeffs', which we embedded in a 32-dimensionsal space with the linear operator 'basis' (a 32x9 matrix).
Since there are only 9 pieces of information in our 32-element signal, we ought to be able to reconstruct the whole thing from nine of its elements.
Let's start by sampling our 32-element signal at 9 (arbitrarily chosen) locations.
To do this, we set up a 9x32 sampling matrix S,
in which each row contains exactly one 1,
in ascending (but otherwise random) columns.
When multiplied by a 32x1 column vector,
the sampling matrix S produces a 9x1 vector selecting 9 of 32 samples.
We create S by deleting 23 randomly chosen rows
from a 32x32 identity matrix:
S = eye(32); for count=1:(32-9), rows = 33-count; r = round(rand*rows); % pick a row at random S = [ S(1:(r-1), :) ; S((r+1):rows, :) ]; % eliminate chosen row end
We can see which 9 rows are selected by S as follows:
>> ( S*(1:32)' )'
ans =
3 9 20 22 23 24 27 28 31
and plot the 9 selected elements of the signal
S_signal = S*signal; plot(1:9,S_signal,'-', 1:9,S_signal,'x');
Now, how to reconstruct the original signal from the 9 samples?
Recall that:
basis is 32x9 random matrix
S is 9x32 sampling matrix
coeffs is 9x1 random coefficient vector
such that signal = basis*coeffs
Since S_signal = S*signal, i.e.
S_signal = S*basis*coeffs,
therefore
inv(S*basis)*S_signal = coeffs
...and of course then we have signal = basis*coeffs
Why not
inv(S)*S_signal = basis*coeffs = signal? Because S is 9x32, not square, and thus does not have an inverse whereas S*basis is a 9x9 square matrix (whose columns contain the sampled elements of the original basis)...
>> recon_coeffs = inv(S*basis) * S_signal;
>> recon_signal = basis * recon_coeffs;
>> [coeffs' ; recon_coeffs']
ans =
Columns 1 through 7
0.6586 0.8636 0.5676 0.9805 0.7918 0.1526 0.8330
0.6586 0.8636 0.5676 0.9805 0.7918 0.1526 0.8330
Columns 8 through 9
0.1919 0.6390
0.1919 0.6390
>> norm(recon_signal-signal)
ans =
6.9783e-14
So it worked -- we got our original signal back from 9 samples.
Could we have done it from 8?
Let's try it. We pick a row of the sampling matrix at random and eliminate it, giving a new 8x32 sampling matrix S2:
>> r = round(rand*9);
>> S2 = [ S(1:(r-1),:) ; S((r+1):9, :) ];
>> ( S2*(1:32)' )'
ans =
3 9 20 22 23 24 28 31
>> ( S*(1:32)' )'
ans =
3 9 20 22 23 24 27 28 31
Now S_signal2 = S2*basis*coeffs -- an 8-element sample of signal -- but S2*basis is not invertible (8x9 instead of 9x9)
What about the pseudo-inverse, which is a general way to undo non-square matrix multiplications, based on SVD.
Recall that if [U S V] = svd(A), for A NxM
so that A = U*S*V'
where U is orthonormal NxN, V is orthonormal MxM, S diagonal NxM
Then the pseudoinverse of A is V*S_ext*U'
where S_ext is an MxN diagonal matrix whose diagonal is
[1./diag(S) zeros(abs(diff(size(S))),1)]
(that is, the inverse of the diagonal of S, as far as it goes,
and then padding with zeros for as many other positions as are needed).
This serves as an "inverse" of A since
V*S_ext*U' * U*S*V' will all cancel appropriately
^^^^^^^^^^ ^^^^^^
pseudoinverse SVD
of A of A
Alternatively and equivalently, we can represent the pseudoinverse
of A as inv(A'*A)*A' , since
inv(A'*A)*A' * A also cancels in the right way.
Both approaches seem to give us a way to reconstruct the signal
from eight samples,
but actually neither will work in this case.
We are trying to get basis*coeffs, or coeffs, back out of the sampled signal, which is S2*basis*coeffs.
We could try for the pseudoinverse of S2, which will involve calculating inv(S2'*S2)*S2', but this is out of the question.
For any sampling matrix S, we will never be able to invert S'*S, since S will always have fewer rows than columns, and so S'*S, whose size is ncolumns(S) x ncolumns(S), will always be full of linear dependencies.
This is just as true for our original 9X32 sampling matrix as for the current 8X32 sampling matrix.
S2*basis will also have fewer rows than columns (8X9 in our example) and so the pseudoinverse construction will also not work here, for the same reason.
We can write down the expressions, and we will get back an answer, but the answer will not be the one we want.
S_signal = S2*signal; M = S2*basis; % M is 8X9 % % pseudoinverse(M) % | % \ / recon_coeffs = inv(M'*M)*M' * S_signal; S_recon = basis*recon_coeffs; S_recon = basis*recon_coeffs;
In fact, we can usually come at least as close by picking a new random combination of the basis vectors -- here we come significantly closer!
norm(signal-S_recon)
ans =
7.1958
norm(signal-basis*rand(9,1))
ans =
5.0711
What about reconstructing the signal from a sample of 10 elements? Since the original signal was constructed from a linear combination of 9 basis vectors, this should be fine.
We need to start with sampling matrix S3, selecting 10 random elements from a 32-element vector:
S3 is a 10x32 matrix
Again, each row contains exactly one 1,
in ascending columns
When S3 is multiplied by a 32x1 column vector,
it produces a 10x1 vector selecting 10 of 32 samples.
S3 = eye(32); for count=1:(32-10), rows = 33-count; r = round(rand*rows); % pick a row at random S3 = [ S3(1:(r-1), :) ; S3((r+1):rows, :) ]; % eliminate chosen row end
The rest is the same as before:
>> ( S3*(1:32)' )'
ans =
1 3 4 5 7 12 24 29 30 31
>> S_signal = S3*signal;
>> M = S3*basis; % M is 10x9
>> % pseudoinverse(M)
>> recon_coeffs = inv(M'*M)*M' * S_signal;
>> S_recon = basis*recon_coeffs;
>> norm(signal-S_recon)
ans =
1.6721e-13
Reconstruction from a subspace of sinusoidal basis vectors
We could do the same thing with sinusoidal (rather than random) basis vectors...
x = (2*pi*(0:31)/32)';
basis = [ ones(32,1) cos(x) sin(x) cos(2*x) sin(2*x) ...
cos(3*x) sin(3*x) cos(4*x) sin(4*x) ];
plot(basis)
signal = basis*coeffs; S_signal = S*signal; plot(1:32,signal,'-',1:32,signal,'o');
>> plot(1:9,S_signal,'-',1:9,S_signal,'o');
>> ( S*(1:32)' )'
ans =
3 9 20 22 23 24 27 28 31
>> recon_coeffs = inv(S*basis) * S_signal; >> recon_signal = basis * recon_coeffs; >> norm(recon_signal-signal) ans = 5.9114e-13
Note that (for this particular basis and this particular sampling matrix) we have a "reconstruction operator"
basis*inv(S*basis)
that will reconstruct a full 32-element vector from 9-element sample.
Just as S is a "downsampling operator" that selects 9 out of 32 vector elements, so basis*inv(S*basis) is an "upsampling operator" that turns a 9-element vector back into a 32-element vector (assuming, of course, that the 9-element vector in questions was sampled from a 32-element vector using S).
We can go back and forth without loss of information as long as we begn with a subspace of 9 or fewer dimensions. If we had started with more than 9 dimensions, then the same downsample-upsample process would yield garbage, instead of giving us back exactly what we started with.
Let's start again with our random linear combination of sinusoidal basis vectors.
signal = basis*coeffs;
What do we predict the magnitude of the dft of 'signal' to be?
DC ((cos(0) -- Fourier coefficient 1) = coeff(1), the cosine and sine terms at frequency 1 = coeff(2) and coeff(3), the cosine and sine terms at frequency 2 = coeff(4) and coeff(5), etc. (with some finesse on the typically annoying scaling issue):
>> coeffs'
ans =
Columns 1 through 7
0.6586 0.8636 0.5676 0.9805 0.7918 0.1526 0.8330
Columns 8 through 9
0.1919 0.6390
>> [2*coeffs(1) abs(coeffs(2:2:9)+coeffs(3:2:9)*i)']
ans =
1.3172 1.0335 1.2603 0.8469 0.6672
>> magspec = abs(fft(signal))/16;
>> magspec(1:5)'
ans =
1.3172 1.0335 1.2603 0.8469 0.6672
>> plot(1:32,magspec,':', 1:32,magspec,'x');
So we get out what we put in, at least once the scale factors are adjusted!
Uniformly sampled signals have some special properties when viewed in the Fourier domain.
Consider what happens when we pick (say) every fourth sample out of a signal. We can view this as two steps: First, set all but every fourth sample to zero; then, "squeeze out" the zeroed samples, leaving only 1/4 of the original sequence. Let's look at the effect of this process one step at a time, working on the same signal constructed out of nine sinusoidal basis vectors.
First, we set 3 of every 4 samples to zero. To do this, multiply pointwise by a signal that has an impulse at every 4th sample and is zero elsewhere.
sampler = zeros(32,1); sampler(1:4:32,1) = ones(8,1); S_signal = sampler.*signal;
We plot the sampled signal:
lplot(S_signal);
Now we plot the spectrum of the sampled signal on the same graph as the spectrum of the original. Notice the difference in amplitude and the induced periodicity -- we now have four copies of the original spectrum:
S_magspec = abs(fft(S_signal))/16; plot(1:32,S_magspec,'xr', 1:32,S_magspec,'og');
Do you see why this is? Having multiplied by the vector sampler in the time domain, we've convolved with the DFT of sampler in the frequency domain. If this isn't clear to you, hold the thought until we get back to it in the general case.
Can we reconstruct the original signal from the partly-zeroed signal?
Well, the original "lives" in a 9-dimensional subspace. The sampled signal contains only eight non-zeroed samples --- 1/4 of 32 = 8 samples were preserved. Since this is informationally equivalent to selecting 8 and throwing the rest away, the previous exercise tells us that we can't reconstruct the original successfully.
Looking at the spectrum, we see the dimensionality problem in a new guise. Our four copies of the original spectrum OVERLAP, and are added together where they do so. This overlap is a form of ALIASING in the frequency domain: two components that should be kept separate are added up, and in general we cannot recover the original operands of an addition if we know only the the result.
In this particular case, because of the complex conjugate relationship of positive- and negative-frequency Fourier coefficients, we can recover the magnitude of the component, but we lose the phase, reflecting the reduction from 9 dimensions to 8.
So as to end on a brighter note, let's try zeroing out only half of the elements of our signal. Since the result preserves 16 of the numbers, the original 9-dimensional signal should be a cinch to reconstruct.
sampler = zeros(32,1); sampler(1:2:32,1) = ones(16,1); S_signal = sampler.*signal; lplot(S_signal);
S_magspec = abs(fft(S_signal))/16; plot(1:32,magspec,':', 1:32,magspec,'x', 1:32,S_magspec,'-',1:32,S_magspec,'o')
Again, by multiplying pointwise by sampler in the time domain, we've created a periodic replication in the frequency domain. This tiime, though, we've only created one extra copy. Now let's reconstitute the original.
Take the DFT -- then correct the scaling factor!
dft = 2*fft(S_signal);
Now let's zero out the replicated copies of the dft coefficients. By doing so we are 'filtering out' unwanted frequency components. This operation is the equivalent of an 'upsample' operator in the time domain, which would be doing the equivalent filtering by convolution.
f_dft = dft; f_dft(10:25) = zeros(1,16);
Now we just reconstitute by inverse DFT:
>> recon = real(ifft(f_dft)); >> norm(signal-recon) ans = 4.6868e-15
What would we get from the inverse DFT without the filtering?
>> plot(1:32,recon1,'-b',1:32,recon1,'ob'); hold on;
>> title('Reconstruction by inverse DFT without filtering')
>> plot(1:32,signal,':r',1:32,signal,'xr');
