How to find spectrum envelope from wav file

How do I measure the energy change of the spectrum envelope 20 - 40 hz from from a wav file and plot it to graph? my recording is in 8000 sampling rate 16bits mono.

I would like to get something like this

Thanks for great help, wish you all merry Christmas and happy new year :)

[Image Analyst]

Daemian, perhaps you overlooked this big hint given in Star's answer "I don’t have your signal, so I can’t test the filter I designed for you with it." HINT, HINT. What do you think you should do now?

Also, please read this.

Answer 1

[Star Strider]

If you have the Signal Processing Toolbox, such a bandpass filter is easy to design. First (to make things easier) convert your 16-bit signed integer signal (that I call ‘x’ here) to double:

xd = double(x);                         % Convert to ‘double’

Fs = 8000;                              % Sampling Frequency
Fn = Fs/2;                              % Nyquist Frequency
Fpb = [20 40]/Fn;                       % Passband
Fsb = [15 50]/Fn;                       % Stopband
[n,Wn] = buttord(Fpb, Fsb, 1, 10);      % Filter Order, With Rp = 1, Rs = 10
[b,a]  = butter(n,Wn);                  % Create Filter Transfer Function
[sos,g] = tf2sos(b,a);                  % Second-Order-Section Implementation

figure(1)
freqz(b,a)                              % Transfer Function Plot

figure(2)                               % Second-Order-Section Plot
freqz(sos)

yd = filtfilt(sos, g, xd);              % Filter ‘xd’ To Get ‘y’

where ‘xd’ is your .wav signal and ‘yd’ is the filtered output. The filtfilt function will filter both channels of your signal at the same time (assuming your mono signal could have two channels, with the same information in both channels), so this code will work regardless of its having 1 or 2 channels. You can then change ‘yd’ to 16-bits with the int16 function if you want to.

Have fun with this!

Merry Christmas (belatedly) and happy New Year to you, too!

[Daemian]

Hi thanks for your reply :) is a great help. but what i'm trying to get the the spectrum envelope like the image below. Do oyou have any idea how do i get it?

[Star Strider]

My pleasure!

I don’t have your signal, so I can’t test the filter I designed for you with it.

I wasn’t sure what you meant by plot (I thought you wanted to plot the transfer function). Again taking ‘yd’ as your filtered signal, this should work:

yd_len = length(yd);                    % Length Of ‘wav’ File Data
Ts = 1/Fs;                              % Sampling Time
Tv = linspace(0, yd_len-1, yd_len)*Ts;  % Time Vector

figure(3)
plot(Tv, yd)                            % Plot Filter Output

You will have to define what you mean by ‘spectrum envelope’ if this doesn’t do what you want. It looks like you want to add your original signal to the output of the filter, in which situation figure(4) here may be what you want:

figure(4)
plot(Tv, yd+xd)                         % Plot Filter Output

This is just a guess.

[26 Dec 2014 - 15:51]

— EDIT —

[28 Dec 2014 - 14:39]

Does this do what you want? If not, what would you want it to do? Please be as specific as possible, and as detailed as necessary.

I have no idea what you mean by ‘energy change of the spectrum envelope’ and can find no online references to anything by that name. If you have a reference, please post the PDF of it, describing the process and what you want to do.

[Daemian]

Hi, i have attach the full pdf for your reference. Is at page 5 section 2.5.1.5 Summary test 1. and the figure i grab from there.

Really appreciate your great help.

Thanks

[Star Strider]

My pleasure.

The problem is that paper gives no information on the signal processing methods it uses to get those plots. It also mentions the physiological tremor frequency band of 8-12 Hz, but does not mention how it interacts with the 20-40 Hz band it mentions elsewhere.

The only reference I can find that has what seems to be a method is behind a paywall (on the Scribd site ‘VOice Stress Analysis’). I do not have access to it.

I have no experience with voice stress analysis. If you provide (and attach here) the ‘VOice Stress Analysis’ paper or another reference that has a detailed description of the methods it uses (in PDF format), I will see if it provides an algorithm I can code in MATLAB.

Without a method or algorithm describing their approach to voice stress analysis, I cannot help. Meanwhile, see if the filter design I posted provides the information you want. If you decide to add a bandpass filter for the 8-12 Hz band, I would design that filter similar to the previous design, only with these passband and stopband changes:

Fpb = [ 8 12]/Fn;                       % Passband
Fsb = [ 5 15]/Fn;                       % Stopband

Everything else is unchanged. Be sure to check the stability of the ‘sos’ implementation of this filter with the freqz function.

[Daemian]

Hi I have attach the voice stress analysis paper see if it help you and figure 3 and 4 and a sample wav file. Really many many thanks :)

[Star Strider]

My pleasure!

Thank you for the paper. It’s 54 pages long and we’re GMT-7 here. I’ll look at it in the morning.

[Daemian]

Sorry about that, Im in Singapore which now 2pm now HAHA!

[Image Analyst]

I guess you chose not to take my hint of attaching your own data. Failing to do that will probably delay your ultimate answer, but that's your choice.

[Daemian]

Sorry i have overlook your comment at the top. Anyway i just find out can attach zip file. I tried to attach a wav alone but unable to. Thanks anyway :)

[Star Strider]

@Image Analyst — Thank you again.

@Daemian — The paper unfortunately does not describe any method for performing the analysis you’re interested in. It goes into detail about the Bayesian decision analysis, but no signal processing. (It refers to Reference [9] for that, apparently.) Please search the literature for a paper that describes in detail the method you want to use. (Section 3.3.2 for instance mentions ‘Praat’.) Then post it here, do your best to code it yourself, and if you need help, we will help you implement it in MATLAB.

Did the filters I designed for you do what you want?

[Daemian]

im not to sure.. after the filters for figure 4 still doesn't look like the Christmas tree waveform . I will take a look at Praat first. Just another question, I found this source http://www.codeproject.com/Articles/7196/ShakyVoice-A-voice-stress-analysis-tool , Is code in C# and is 10 years ago code. do you have any idea is it possible to implement this in matlab? I tried to implement but im not too sure how he get the value.

[Star Strider]

It is possible to implement C and FORTRAN and other compiled code as MATLAB mex files, but I’ve never done it. (I never needed to. I did a lot of FORTRAN programming, but never implemented any as mex functions.)

I am certain that MATLAB can do everything you want. You simply need to identify the algorithm you want to implement. Just now, those appear to be proprietary, so I am not certain how you will get them, but there may be a reference somewhere that describes them. I suggest you do a PubMed search. One of those references may provide you with the information you need.

[Daemian]

Thanks a lot! anyway HAPPY NEW YEAR!

[Star Strider]

My pleasure!

Happy New Year to you, too!

[Daemian]

Hey Star Strider, I found this pdf talking about spectral envelopes I not too sure is it helpful, can you help me take a look?

Thanks!

Answer 2

[Youssef Khmou]

Hilbert Transformation is used to obtain the envelope of signal, here is an example :

 The envelope is decreasing exponential, 

Fs=80;
F=10;
t=0:1/Fs:4-1/Fs;
x=exp(-t).*real(exp(j*2*pi*F*t));
figure; plot(t,x);
Y=abs(hilbert(x));
hold on;
plot(t,Y,'r');

 fx=fftshift(abs(fft(x))); fx=fx(floor(end/2:end));
 fY=fftshift(abs(fft(Y))); fY=fY(floor(end/2:end));
 figure; plot(fx); hold on
 plot(fY,'r')

[Image Analyst]

Cool - I didn't know hilbert() did that. Does it always get the envelope no matter what x looks like (like how fast it oscillates)? Why does it get a little squirrely around t=4?

[Youssef Khmou]

According to theory yes, oscillation in borders can be interpreted as Gibbs effect.