fidelity_error.m

来自「fastest algorithm to find EMD.」· M 代码 · 共 63 行

function fidelity_error(freq,Ny_multiple);
% fidelity_error(freq,Ny_multiple);
% 
% This script demonstrates how sampling faster than the Nyquist
% frequency forces amplitude modulation for all sample frequencies
% that are not integer multiples of the signal's frequency.  We consider
% this behavior to be a loss-of-fidelity error.  It is not aliasing because 
% the frequency content is not changed, just the amplitude. You will 
% see that the error is reduced as the sample frequency increases.
% 
% Inputs:
% freq          = frequency to test in Hz
% Ny_multiple   = how many times faster than Nyquist to model to
%
% example: to model to 5 times Nyquist for a 1 Hz tone,
%   fidelity_error(1,5);
% 
% It may seem nonsensical to even try to sample a given pure tone at a rate
% that is anything but an integer multiple of that frequency, but this is 
% exactly what happens with seismic pulses of infinite bandwidth.  There 
% exists no pure tone in a seismic pulse with which we can decide the "proper"
% sample frequency.  Therefore, any loss of information is going to produce 
% fidelity errors when using seismic pulses, and the only way to
% minimize the error is to sample as fast as possible
%
% author: Bradley Matthew Battista
%   University of South Carolina
%   Department of Geological Sciences
%   701 Sumter Street, EWS 617
%   Columbia, SC. 29208
%
% COPYRIGHT: see the associated COPYRIGHT.txt file, and also
% http://software.seg.org/disclaimer2.txt
% This source code may be found online at:
% http://software.seg.org/2007/0003
%

% use 20*Nyquist for comparison sample frequency
tt = 0:(1/(40*freq)):1;
xx = sin(2*pi*freq*tt);
plot(tt,xx,'-b');hold on
h = plot(1,1,'.');

% loop thru sample frequencies Nyquist thru Ny_multiple
% and compare the time series to the control signal above
for n = (2*freq*[1:.01:Ny_multiple])
    delete(h);
    dt = 1/n;              % determine sample interval
    t  = 0:dt:1;           % generate time samples
    x  = sin(2*pi*freq*t); % calculate amplitude for each sample
    h = plot(t,x,'.-r');    % plot the time-series
    title(['Nyquist \times ',num2str(n/(2*freq))],'fontsize',14); 
    set(gca,'ylim',[-1 1]);drawnow;pause(.1);
    
    % for a sine, the number of peaks and troughs
    % should equal 2*freq.  There should be 2*freq
    % samples with |amp|=1 if no loss of fidelity is occurring
    % Lets display a message when this occurs
    if length(find(abs(x) == 1)) == 2*freq
        disp(['Nyquist x ',num2str(n/(2*freq)),', no error']);
    end
end