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 frequencytt = 0:(1/(40*freq)):1;xx = sin(2*pi*freq*tt);plot(tt,xx,'-b');hold onh = plot(1,1,'.');% loop thru sample frequencies Nyquist thru Ny_multiple% and compare the time series to the control signal abovefor 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']);    endend