multistage sample-凯发k8网页登录

setup

define some parameters to be used throughout the example.

framesize = 64;
infs      = 96e3;

generating the 96 khz signal

generate the chirp signal using dsp.chirp as follows:

source = dsp.chirp(initialfrequency=0,targetfrequency=48e3, ...
    sweeptime=8,targettime=8,samplerate=infs, ...
    samplesperframe=framesize,type="quadratic");

create spectrum analyzers

create two spectrum analyzers. these will be used to visualize the frequency content of the original signal as well as that of the signals converted to 44.1 khz.

spectrumanalyzer44p1 = spectrumanalyzer( ...
    samplerate=44100, ...
    method="welch", ...
    averagingmethod="exponential", ...
    forgettingfactor=1e-7, ...
    viewtype="spectrum-and-spectrogram", ...
    timespansource="property",timespan=8, ...
    window="kaiser",sidelobeattenuation=220, ...
    ylimits=[-250, 50],colorlimits=[-150, 20], ...
    plotastwosidedspectrum=false);
spectrumanalyzer96 = spectrumanalyzer( ...
    samplerate=96000, ...
    method="welch", ...
    averagingmethod="exponential", ...
    forgettingfactor=1e-7, ...
    viewtype="spectrum-and-spectrogram", ...
    timespansource="property",timespan=8, ...
    window="kaiser",sidelobeattenuation=220, ...
    ylimits=[-250, 50],colorlimits=[-150, 20], ...
    plotastwosidedspectrum=false);

spectrum of original signal sampled at 96 khz

the loop below plots the spectrogram and power spectrum of the original 96 khz signal. the chirp signal starts at 0 and sweeps to 48 khz over a simulated time of 8 seconds.

nframes = 8*infs/framesize;
for k = 1:nframes
    sig96 = source();          % source
    spectrumanalyzer96(sig96); % spectrogram
end
release(source)
release(spectrumanalyzer96)

setting up the sample rate converter

in order to convert the signal, dsp.samplerateconverter is used. a first attempt sets the bandwidth of interest to 40 khz, i.e. to cover the range [-20 khz, 20 khz]. this is the usually accepted range that is audible to humans. the stopband attenuation for the filters to be used to remove spectral replicas and aliased replicas is left at the default value of 80 db.

bw40 = 40e3;
outfs = 44.1e3;
src40khz80db = dsp.samplerateconverter(bandwidth=bw40, ...
    inputsamplerate=infs,outputsamplerate=outfs);

analysis of the filters involved in the conversion

use info to get information on the filters that are designed to perform the conversion. this reveals that the conversion will be performed in two steps. the first step involves a decimation by two filter which converts the signal from 96 khz to 48 khz. the second step involves an fir rate converter that interpolates by 147 and decimates by 160. this results in the 44.1 khz required. the freqz command can be used to visualize the combined frequency response of the two stages involved. zooming in reveals that the passband extends up to 20 khz as specified and that the passband ripple is in the milli-db range (less than 0.003 db).

info(src40khz80db)
[h80db,f] = freqz(src40khz80db,0:10:25e3);
plot(f,20*log10(abs(h80db)/norm(h80db,inf)))
xlabel("frequency (hz)")
ylabel("magnitude (db)")
axis([0 25e3 -140 5])

ans =
    'overall interpolation factor    : 147
     overall decimation factor       : 320
     number of filters               : 2
     multiplications per input sample: 42.334375
     number of coefficients          : 8618
     filters:                         
        filter 1:
        dsp.firdecimator     - decimation factor   : 2 
        filter 2:
        dsp.firrateconverter - interpolation factor: 147
                             - decimation factor   : 160 
     '

asynchronous buffer

the sample rate conversion from 96 khz to 44.1 khz produces 147 samples for every 320 input samples. because the chirp signal is generated with frames of 64 samples, an asynchronous buffer is needed. the chirp signal is written 64 samples at a time, and whenever there are enough samples buffered, 320 of them are read and fed to the sample rate converter.

buff = dsp.asyncbuffer;

main processing loop

the loop below performs the sample rate conversion in streaming fashion. the computation is fast enough to operate in real time if need be.

the spectrogram and power spectrum of the converted signal are plotted. the extra lines in the spectrogram correspond to spectral aliases/images remaining after filtering. the replicas are attenuated by better than 80 db as can be verified with the power spectrum plot.

srcframesize = 320;
for k = 1:nframes
    sig96 = source();       % generate chirp
    write(buff,sig96);      % buffer data
    if buff.numunreadsamples >= srcframesize
        sig96buffered = read(buff,srcframesize);
        sig44p1 = src40khz80db(sig96buffered); % convert sample-rate
        spectrumanalyzer44p1(sig44p1);   % view spectrum of converted signal
    end
end
release(source)
release(spectrumanalyzer44p1)
release(buff)

a more precise sample rate converter

in order to improve the sample rate converter quality, two changes can be made. first, the bandwidth can be extended from 40 khz to 43.5 khz. this in turn requires filters with a sharper transition. second, the stopband attenuation can be increased from 80 db to 160 db. both these changes come at the expense of more filter coefficients over all as well as more multiplications per input sample.

bw43p5 = 43.5e3;
src43p5khz160db = dsp.samplerateconverter(bandwidth=bw43p5, ...
    inputsamplerate=infs,outputsamplerate=outfs, ...
    stopbandattenuation=160);

analysis of the filters involved in the conversion

the previous sample rate converter involved 8618 filter coefficients and a computational cost of 42.3 multiplications per input sample. by increasing the bandwidth and stopband attenuation, the cost increases substantially to 123896 filter coefficients and 440.34 multiplications per input sample. the frequency response reveals a much sharper filter transition as well as larger stopband attenuation. moreover, the passband ripple is now in the micro-db scale. note: this implementation involves the design of very long filters which takes several minutes to complete. however, this is a one time cost which happens offline (before the actual sample rate conversion).

info(src43p5khz160db)
[h160db,f] = freqz(src43p5khz160db,0:10:25e3);
plot(f,20*log10(abs(h160db)/norm(h160db,inf)));
xlabel("frequency (hz)")
ylabel("magnitude (db)")
axis([0 25e3 -250 5])

ans =
    'overall interpolation factor    : 147
     overall decimation factor       : 320
     number of filters               : 2
     multiplications per input sample: 440.340625
     number of coefficients          : 123896
     filters:                         
        filter 1:
        dsp.firdecimator     - decimation factor   : 2 
        filter 2:
        dsp.firrateconverter - interpolation factor: 147
                             - decimation factor   : 160 
     '

main processing loop

the processing is repeated with the more precise sample rate converter.

once again the spectrogram and power spectrum of the converted signal are plotted. notice that the imaging/aliasing is attenuated enough that they are not visible in the spectrogram. the power spectrum shows spectral aliases attenuated by more than 160 db (the peak is at about 20 db).

for k = 1:nframes
    sig96 = source();              % generate chirp
    over = write(buff,sig96);      % buffer data
    if buff.numunreadsamples >= srcframesize
        [sig96buffered,under] = read(buff,srcframesize);
        sig44p1 = src43p5khz160db(sig96buffered); % convert sample-rate
        spectrumanalyzer44p1(sig44p1);   % view spectrum of converted signal
    end
end
release(source)
release(spectrumanalyzer44p1)
release(buff)