sound sample too short to see higher frequencies?
2 views (last 30 days)
Show older comments
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')
0 Comments
Accepted Answer
Dr. Seis
on 5 Oct 2011
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.
2 Comments
Dr. Seis
on 5 Oct 2011
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.
More Answers (1)
Rick Rosson
on 5 Oct 2011
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
2 Comments
See Also
Categories
Find more on Pulsed Waveforms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!