www.gusucode.com > signal 案例源码程序 matlab代码 > signal/ReconstructingMissingDataExample.m

    %% Reconstructing Missing Data
% This example shows how to reconstruct missing data via interpolation,
% anti-aliasing filtering, and autoregressive modeling.

%% Introduction
% With the advent of cheap data acquisition hardware, you often have access
% to signals that are rapidly sampled at regular intervals.  This allows
% you to gain a fine approximation to the underlying signal.  But what
% happens when the data you are measuring are coarsely sampled or otherwise
% missing significant portions?  How do you infer the values of the signals
% at points in between the samples that you know?

%% Linear Interpolation
% Linear interpolation is by far the most common method of inferring values
% between sampled points.  By default, when you plot a vector in MATLAB,
% you see the points connected by straight lines.  You need to sample a
% signal at very fine detail in order to approximate the true signal.
%
% In this example, a sinusoid is sampled with both fine and coarse
% resolution.  When plotted on a graph, the finely sampled sinusoid very
% closely resembles what the true continuous sinusoid would look like.
% Thus, you can use it as a model of the "true signal." In the plot
% below, the samples of a coarsely sampled signal are shown as circles
% connected by straight lines.

tTrueSignal = 0:0.01:20;
xTrueSignal = sin(2*pi*2*tTrueSignal/7);

tSampled = 0:20;
xSampled = sin(2*pi*2*tSampled/7);

plot(tTrueSignal,xTrueSignal,'-', ...
    tSampled,xSampled,'o-')
legend('true signal','samples')

%%
% It is straightforward to recover intermediate samples in the same way
% that |plot| performs interpolation.  This can be accomplished with the
% linear method of the |interp1| function.

tResampled = 0:0.1:20;
xLinear = interp1(tSampled,xSampled,tResampled,'linear');
plot(tTrueSignal,xTrueSignal,'-', ...
     tSampled, xSampled, 'o-', ...
     tResampled,xLinear,'.-')
legend('true signal','samples','interp1 (linear)')

%%
% The problem with linear interpolation is that the result is not very
% smooth.  Other interpolation methods can produce smoother approximations.

%% Spline Interpolation
% Many physical signals are like sinusoids in that they are continuous and
% have continuous derivatives.  You can reconstruct such signals by using
% cubic spline interpolation, which ensures that the first and second
% derivatives of the interpolated signal are continuous at every data point:

xSpline = interp1(tSampled,xSampled,tResampled,'spline');
plot(tTrueSignal,xTrueSignal,'-', ...
     tSampled, xSampled,'o', ...
     tResampled,xLinear,'.-', ...
     tResampled,xSpline,'.-')
legend('true signal','samples','interp1 (linear)','interp1 (spline)')

%%
% Cubic splines are particularly effective when interpolating signals that
% consist of sinusoids.  However, there are other techniques that can be
% used to gain greater fidelity to physical signals which have continuous
% derivatives up to a very high order.

%% Resampling with Antialiasing Filters
% The |resample| function in the Signal Processing Toolbox provides another
% technique to fill in missing data. |resample| can reconstruct sinusoidal
% components of low frequency with very low distortion.

xResample = resample(xSampled, 10, 1);
tResample = 0.1*((1:numel(xResample))-1);

plot(tTrueSignal,xTrueSignal,'-', ...
     tResampled,xSpline, '.', ...
     tResample, xResample,'.')
legend('true signal','interp1 (spline)','resample')

%%
% Like the other methods, |resample| has some difficulty reconstructing the
% endpoints.  On the other hand, the central portion of the reconstructed
% signal agrees very well with the true signal.
xlim([6 10])


%% Resampling with Missing Samples
% |resample| can accommodate nonuniformly sampled signals.  The technique
% works best when the signal is sampled at a high rate.
%
% In the following example, we create a slowly moving sinusoid, remove
% a sample, and zoom into the vicinity of the missing sample.

tTrueSignal = 0:.1:20;
xTrueSignal = sin(2*pi*2*tTrueSignal/15);

Tx = 0:20;

Tmissing = Tx(10);
Tx(10) = [];

x = sin(2*pi*2*Tx/15);
Xmissing = sin(2*pi*2*Tmissing/15);

[y, Ty] = resample(x,Tx,10,'spline');

plot(tTrueSignal, xTrueSignal, '-', ...
     Tmissing,Xmissing,'x ', ...
     Tx,x,'o ', ...
     Ty,y,'. ')
legend('true signal','missing sample','remaining samples','resample with ''spline''')
ylim([-1.2 1.2])
xlim([6 14])
%%
% The reconstructed sinusoid tracks the shape of the true signal reasonably
% well, with only a slight error in the vicinity of the missing sample.
%
% However, |resample| does not work well when there is a large gap in the
% signal.  For example, consider a dampened sinusoid whose middle portion
% is missing:

tTrueSignal = (0:199)/199;
xTrueSignal = exp(-0.5*tTrueSignal).*sin(2*pi*5*tTrueSignal);

tMissing = tTrueSignal;
xMissing = xTrueSignal;
tMissing(50:140) = [];
xMissing(50:140) = [];

[y,Ty] = resample(xMissing, tMissing, 200, 'spline');

plot(tTrueSignal,xTrueSignal,'-', ...
     tMissing,xMissing,'o',...
     Ty,y,'.')
legend('true signal','samples','resample with ''spline''')

%%
% Here |resample| ensures that the reconstructed signal is continuous and
% has continuous derivatives in the vicinity of the missing points.
% However, it cannot adequately reconstruct the missing portion.  

%% Reconstructing Large Gaps
% As can be seen above, filtering and cubic interpolation alone might not
% be sufficient to deal with large gaps.  However, with certain kinds of
% sampled signals, such as those that arise when observing oscillating
% phenomena, you can often infer the values of missing samples based upon
% the data immediately preceding or following the gap. 
%
% The |fillgaps| function can replace missing samples (specified by |NaN|)
% in an otherwise uniformly sampled signal by fitting an autoregressive
% model to the samples surrounding a gap and extrapolating into the gap
% from both directions.

tTrueSignal = (0:199)/199;
xTrueSignal = exp(-.5*tTrueSignal).*sin(2*pi*5*tTrueSignal);
gapSignal = xTrueSignal;
gapSignal(50:140) = NaN;
y = fillgaps(gapSignal);

plot(tTrueSignal,xTrueSignal,'-', ...
     tTrueSignal,gapSignal,'o', ...
     tTrueSignal,y,'.')

legend('true signal','samples','reconstructed signal')

%%
% The technique works because autoregressive signals have information that
% is spread out over many samples.  Only a few samples within any segment
% are needed to completely reconstruct the full signal.

%%
% This type of reconstruction can be adapted to estimate missing samples of
% more complicated signals.  Consider a sampled audio signal of a plucked
% guitar string after removing six hundred samples immediately after the
% pluck:

[y,fs] = audioread('guitartune.wav');
x = y(1:3500);
x(2000:2600) = NaN;
y2 = fillgaps(x);
plot(1:3500, y(1:3500), '-', ...
     1:3500, x, '.', ...
     1:3500, y2, '-')
legend('original signal','samples','reconstructed signal',...
       'Location','best')

%% Reconstructing Gaps with Localized Estimation
% It is fairly straightforward to fill data within a gap when it is known
% that the signal in the vicinity of the gap can be modeled with a single
% autoregressive process.  You can mitigate problems when a signal consists
% of a non-constant autoregressive process by restricting the area over
% which the model parameters are computed.  This is useful when you are
% trying to fill a gap within the "ringing" period of one resonance that
% comes immediately before or after another, stronger, resonance.

x = y(350001:370000);
x(6000:6950) = NaN;
y2 = fillgaps(x);
y3 = fillgaps(x,1500);
plot(1:20000, y(350001:370000), '-', ...
     1:20000, x, '.', ...
     1:20000, y2, '-', ...
     1:20000, y3, '-')
xlim([2200 10200])
legend('original signal','samples','fillgaps (all)','fillgaps (localized)',...
       'Location','best')

%%
% In the plot above, the waveform is missing a section just before a large
% resonance.  As before, |fillgaps| is used to extrapolate into the gap
% region using all available data.  A second call to |fillgaps| uses only
% 1500 samples on either side of the gap to perform the modeling.  This
% mitigates the effect of the subsequent guitar pluck after sample 7500.

%% Summary
% You have seen several ways to reconstruct missing data from its
% neighboring sample values using interpolation, resampling and
% autoregressive modeling.  
% 
% Interpolation and resampling work for slowly varying signals.  Resampling
% with antialiasing filters often does a better job at reconstructing
% signals that consist of low-frequency components.  For reconstructing
% large gaps in oscillating signals, autoregressive modeling in the
% vicinity of the gap can be particularly effective.
     
%% Further reading
% For more information on uniform and non-uniform resampling, see
%
% * <http://www.mathworks.com/help/signal/examples/resampling-uniformly-sampled-signals.html
% Resampling Uniformly Sampled Signals>
% * <http://www.mathworks.com/help/signal/examples/resampling-nonuniformly-sampled-signals.html
% Resampling Nonuniformly Sampled Signals>