sound sample too short to see higher frequencies?

High Frequency FFT

Hi,

I have a Wav file with a length of approximately 292kB at 88bt/sec, for an old lister engine. I want to change the variable FreqMax to look at frequencies up to the rated RPM (around 1640), but am not able to do this because I don’t have enough data points from this file. (I understand that for 88kbps, my nyquist frequency will not be high enough to get too information, but I want to anyway.) I have three ideas about how to accomplish this but I am unsure how to implement them:

1.  Reduce the frequency resolution

2.  Drop the lower frequencies, to shift the range toward the higher frequencies

3.  Zero-pad the original sample so that I have more data points.

Are there better solutions?

I am curious how (1) or (2) can be accomplished, I am unsure how to do the scaling and cannot find good resources on this, and also if there is some better approach to looking at frequencies with this signal.

Thanks, Chris

%lister engine sound source
%http://www.oldengineshed.com/diesel.html
%3 3/16" Bore X 3 1/2" Stroke - Rated H.P. at 1650 R.P.M.

%FULL SOUND FILE
[x1 fs1 nbits1 opts1]=wavread('listersl.wav');
t1=(0:length(x1)-1)/fs1;
subplot(5,1,1);
plot(t1,x1)

legend('Full Dataset Waveform');
xlabel('t (s)');
ylabel('Amplitude');

%ANALYSIS OF SAMPLE WINDOW
%file bit rate
bitrate=88000/8; %bit/s (8bit/byte)
SampleStart=10; %what second in the file to start reading
SampleEnd=18;
WavSampleLength=[SampleStart*bitrate+1 SampleEnd*bitrate];

[x fs]=wavread('listersl.wav', WavSampleLength);

t=(SampleStart:1/fs:SampleEnd);

dif=length(t)-length(x);

t=(SampleStart:1/fs:SampleEnd-1/fs*dif);

subplot(5,1,2);
plot(t,x)

legend('Sample Waveform');
xlabel('t (s)');
ylabel('Amplitude');

xwin=x.*hamming(length(x));

subplot(5,1,3)
plot(t,xwin)

legend('Window');
xlabel('t (s)');
ylabel('Amplitude');

Y=fft(xwin);
FreqMax=10000;
hz=FreqMax*length(Y)/fs;
f=(0:hz)*fs/length(Y);
subplot(5,1,4);
plot(f,20*log10(abs(Y(1:length(f)))+eps));
legend('Spectrum');
xlabel('Freq (Hz)')
ylabel('Magnitude (dB)')

C=fft(log(abs(Y)+eps));
f1=fs/1000; %max for plot, fs/1000~1 milliseconds
f2=fs/1; %min for plot, fs/50~20 milliseconds

q=(f1:f2)/fs;
subplot(5,1,5);
plot(q,abs(C(f1:f2)));
legend('Cepstrum')
xlabel('Quefrequency (s)')
ylabel('Amplitude')

Answer 1

[Dr. Seis]

The Nyquist frequency is, by definition, the highest frequency you can resolve before you start aliasing the signal. This means the amplitude you derive for any frequency beyond the Nyquist will likely have little meaning. You find the Nyquist by:

sample_rate = 20; % units = samples per second;

nyquist_freq = sample_rate / 2; % = 10 Hz

Notice this equation has nothing to do with how long the timeseries is. The length of timeseries along with the sample rate controls the frequency increment (the smallest resolvable frequency) like this:

sample_rate = 20;
num_samples = 1024;
time_increment = 1 / sample_rate;
freq_increment = sample_rate / num_samples;

The fact is, there is no way to recover or resolve higher frequencies unless you have a higher sample rate. Period. If you can't increase the sample rate of your measuring device and you want bogus results, then you can interpolate between the points in your time series to increase the sample rate. You need to interpolate to a sample rate that pushes the Nyquist frequency beyond the highest frequency you will be resolving bogus-ly.

[Chris]

Thanks Elige,

That part of it I understand. I was hoping to ask more about sample length rather than rate - so in a situation where the sample rate is good enough, but the length of the sample is very short, it seems like it should be readily possible to see higher frequency responses but I'm not sure how to get at them from the examples I've looked at.

Hypothetically, if this file were shorter in length but with the same kbps, could I still detect the frequencies I currently detect in my FFT(x) vs Frequency, and if so, how? Currently, I would have to lower FreqMax so that f=(0:hz)*fs/length(Y) and Y=fft(xwin) could still have the same length.

[Dr. Seis]

So, if I am gathering the info in all the comments here correctly, the Nyquist frequency of your data is 5512.5 Hz and the frequency you are trying to resolve is 27.33 Hz (1640 RPM). So, if you don't have enough time samples, then your frequency increment will be a lot higher (per the example above) than you want. If this is the case, then padding the timeseries with zeros will increase your "num_samples" and decrease your "freq_increment" above. To get down to 27.33 Hz, you will need at least 404 samples (= 11025 samples per second / 27.33 Hz) in the time domain (preferably a power of two like 2^9 = 512). I would pad evenly around both sides of the data you do have, instead of only one side. It may also be necessary to create a taper where the zeros meet your data, so that there isn't a sharp response in the frequency domain.

Answer 2

[Rick Rosson]

Hi Chris,

What is the value of the sampling rate fs (in samples per second) as returned in the first line of code?

   [x fs] = waveread(...);

The bit rate (in bits per second) is not relevant to this discussion.

Thanks!

Rick

[Chris]

fs =11025.

[Rick Rosson]

Okay, then you can represent frequencies up to fs/2 = 5512.5 hertz.