5
$\begingroup$

Im reading a book about DSP and there is an example about investigating the sound that travel through the ocean. The sound was recorded as the time-domain signal and they tried to analyze it in frequency domain.

The input signal was broke up into 256 point segments. Each of these segments is multiplied by the Hamming window, then run through a 256 point DFT. The resulting frequency spectra are then averaged to form a single 129 point frequency spectrum. Please take a look at the figure

enter image description here

I don't understand:

  1. How they break up the input signal into 256 segments?

  2. How each segment was multiplied with the Hamming window?

  3. How they calculate 256 DFT points?

  4. How the frequency spectra was averaged?

$\endgroup$

migrated from stackoverflow.com Aug 27 '12 at 16:20

This question came from our site for professional and enthusiast programmers.

9
$\begingroup$

How they break up the input signal into 256 segments?

This is usually done by a sliding a window of a given width (here 256). Depending on the step-size this results in overlapping (if smaller than 256) or non-overlapping (if equal or larger than 256) data segments.

How each segment was multiplied with the Hamming window?

Usually, by an element-wise multiplication.

How they calculate 256 DFT points?

The number of DFT points is, if I am not mistaken, determined by the length of the input signal.

How the frequency spectra was averaged?

All resulting spectra were summed up and divided by the number of spectra.

$\endgroup$
  • $\begingroup$ Thanks Muster, but to be honest, I have not yet understood what you've explained (sorry.) I'm a student and brandnew to DSP :( so could you or anybody plz help to explain to me in more detail ? or could you tell me the material with which I can read to understand? $\endgroup$ – Mai Aug 27 '12 at 14:41
  • $\begingroup$ @Mai With all due respect, if you need more detail than this you should probably buy a signal processing book and study it or take a class. $\endgroup$ – Jim Clay Aug 27 '12 at 18:23
  • $\begingroup$ I am learning myself by eading this book "The Scientist and engineer's Guide to Digital Signal Processing" and stuck here :( OK, I will try to figure it out. Thank you :) $\endgroup$ – Mai Aug 27 '12 at 20:37
  • 1
    $\begingroup$ @Mai My apologies then. If you install Octave (a freeware program that is a lot like Matlab) on your computer you should be able to run the code that Christian has posted. It is doing the steps that you are asking about. Learn how Octave works, run the code and examine the results and how it works, and you should get the answers that you are looking for. $\endgroup$ – Jim Clay Aug 28 '12 at 0:57
3
$\begingroup$

The code for this could look something like this (assuming you're only interested in magnitude informations of the spectrum):

clear all;
clc;

input_signal = sin(2*pi*.2*[1:256*1024]).'; %sine at frequency .2 (n=51.2)'

for idx = 1:length(input_signal)/256
     %take a slice from the input vector
    slice = input_signal(1+(idx-1)*256:idx*256);
    %multiply it with the window and transform it into frequency domain
    spectrum = fft(slice.*hamming(256));
    %get the spectrum magnitude at each of the 256 frequency points and store it
    mag_spectrum(:,idx) = abs(spectrum).^2;
end

%for each frequency, take the mean value of the magnitude (mean along rows, thus the .')
mean_spectrum = mean(mag_spectrum.');

After that, the vector mean_spectrum contains positive and negative frequencies (positive frequencies at indexes 2..128, negative frequencies (-fs/2 first!) at indexes 129..256, DC at index 1). Thus, I assume the output vector of length 128 averages positive and negative frequencies. You could obtain this by:

output_vec(1) = mag_spectrum(1);
output_vec(2:128) = 1/2*(mag_spectrum(2:128)+mag_spectrum(256:-1:130));
output_vec(129) = mag_spectrum(129);
plot(output_vec);

This results in a peak in the spectrum at approx. n = 52, as expected. If you want to track the frequency components of the signal over time, you might use a Short-Time-Fourier-Transformation instead.

$\endgroup$
  • $\begingroup$ Thanks Christian, I got it :) I have a time-based signal sampled at 6MHz. If I want to have frequency resolution of about 3Hz then I have to use 2^21 FFT points. Am I correct? Could you please take a look at the following code and help to check if it is correct: L=2^21 for index=1:length(input)/L slice = input(1+(index-1)*L:indexL); abcd = han(L) . slice'; NFFT = 2^nextpow2(L); spectrum = fft(abcd, NFFT)/L; mag_spectrum(:, index)=abs(spectrum.^2); end $\endgroup$ – Mai Aug 28 '12 at 12:24
  • $\begingroup$ oh, sorry. Could any Mod help to delete the first comment plz? Hope you can see the code :( $\endgroup$ – Mai Aug 28 '12 at 12:26
  • 1
    $\begingroup$ @Mai Yes, 2^21 points will get you 3 Hz resolution at 6 Msamples/s. If you want people to review your code, you should probably ask a question at Stackoverflow. $\endgroup$ – Jim Clay Aug 28 '12 at 12:59
  • $\begingroup$ @JimClay stackoverflow? Are you suggesting DSP'ish code snips should be asked on stack overflow and not dsp.se? I would think a separate question, specific to the code, on dsp.se linking to the original question would be a good idea. $\endgroup$ – Christopher Felton Aug 29 '12 at 12:20
  • 1
    $\begingroup$ @JimClay I would disagree. I would think Signal processing code and Signal processing topics should be posted here. Obviously, some grey area here, if someone wanted to implement their own PSD and have code issues, I would think dsp.se is appropriate. If someone had a generic question, how to write a function in Matlab, overflow might be more appropriate? I think it would still be ok on dsp.se as well. $\endgroup$ – Christopher Felton Aug 29 '12 at 14:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.