matlab rectify,matlab 声音分析[转]

Lecture 9: Spectral Analysis

Objectives

By the end of the session you should:

qknow

how to build and perform a filterbank analysis of a signal

qknow

how to use the discrete fourier transform to convert a waveform to

a spectrum and vice versa.

qknow

how to divide a signal into overlapping windows

qknow

how to calculate and display a spectrogram

Outline

1. Filterbank

analysis

The most flexible way to

perform spectral analysis is to use a bank of bandpass

filters. A filterbank can be designed to provide

a spectral analysis with any degree of frequency resolution (wide

or narrow), even with non-linear filter spacing and

bandwidths. A disadvantage of filterbanks is that

they almost always take more calculation and processing time than

discrete Fourier analysis using the FFT (see below).

To use a filter bank for

analysis we need one band-pass filter per channel to do the

filtering, a means to perform rectification, and a low-pass filter

to smooth the energies. In this example, we build

a 19-channel filterbank using bandwidths that are modelled on human

auditory bandwidths. We rectify and smooth the

filtered energies and convert to a decibel scale.

function

e=auditoryfbank(x,fs)

% AUDITORYFBANK performs an

auditory filterbank analysis on a signal

% E=AUDITORYFBANK(X,FS)

passes signal X, sampled at FS samples/sec

% through a 19-channel

auditory filterbank, smoothing and downsampling

% energies to 100

frames per second. Output is in decibels.

%

% get output energy sample

points

x=reshape(x,[length(x)

1]); % force to be a column

eidx=fs/100:fs/100:length(x); % sampling indices (100 frame/sec)

ne=length(eidx); % number of output energy frames

% these are the filter

cut-offs

cuts=[180 300 420 540 660 780

900 1075 1225 1375 1525 1700 1900 2100 2300 2550 2850 3150 3475

4000];

nf=length(cuts)-1; % number of filters

% initialise the output energy

table

e=zeros(ne,nf);

% build a smoothing

filter

[sb,sa]=butter(4,2*50/fs); % low-pass at half output frame rate

% for each channel in

turn

for i=1:nf

fprintf('Calculating channel %d

(%d-%dHz)\n',i,cuts(i),cuts(i+1));

%

build the band-pass filter

[b,a]=butter(4,[2*cuts(i)/fs 2*cuts(i+1)/fs]);

%

filter the signal

y=filter(b,a,x);

%

rectify and smooth

sy=filter(sb,sa,abs(y));

%

sample and convert to decibels

sy=max(sy,eps); % force to be > 0

e(:,i)=100+20*log10(sy(eidx)); % downsample and convert to dB

end

The input to this function is a

waveform,x, sampled

atfssamples/second. The

output is a table of energies,e, in which each row represents 10ms of the

signal, and each column the energy in decibels found in a band-pass

filtered frequency region at that time.

In this example, we build a

chirp and pass it through the filterbank, then plot the output

energies in each channel against time:

% build a chirp from 500 to

2000Hz

t=0:1/10000:1;

f=linspace(500,2000,length(t));

y=sin(2*pi*f.*t);

soundsc(y,10000);

% put through filter

bank

e=auditoryfbank(y,10000);

% scale to 0..1, and

transpose

emin=min(min(e));

emax=max(max(e));

te=((e-emin)/(emax-emin))'

% offset each channel by

channel number for plot

for i=1:size(te,1)

te(i,:) = te(i,:)+i;

end

plot(1:size(te,2),te,'k-');

ylabel('Filter

channel');

xlabel('Frame

number');

2. Spectral analysis using Fourier transform

The discrete-time discrete-frequency version of the Fourier

transform (DFT) converts an array of N sample amplitudes to an

array of N complex harmonic amplitudes. If the

sampling rate is Fs, the N input samples are 1/Fs seconds apart,

and the output harmonic frequencies are Fs/N hertz

apart. That is the N output amplitudes are evenly

spaced at frequencies between 0 and (N-1).Fs/N hertz.

To compute the DFT in MATLAB, we use the

function fft(x,n). This function takes a waveform x and the number of samples

n. When n is less than the length of x, then x is

truncated; when n is longer than the length of x, then x is padded

with zeros. The output is an array of complex

amplitudes of length n. You can obtain the

magnitude of each spectral component

with abs(), and its

phase with angle() (result

in radians).

In this example, we generate a complex periodic signal, then

calculate and display the magnitude spectrum and the phase

spectrum:

% generate a complex signal

fs =

10000; %

sampling rate

t =

0:1/fs:0.1; % sampling instants

x = sin(2*pi*1500*t) +

sin(2*pi*4000*t); % two sinusoids

%

% perform 1000-point transform

y =

fft(x,1000); % y contains 1000 complex amplitudes

y =

y(1:500); % just look at first

half

m =

abs(y); % m = magnitude of sinusoids

p =

unwrap(angle(y)); % p = phase of sinusoids, unwrap()

% copes with 360 degree jumps

% plot spectrum 0..fs/2

f =

(0:499)*fs/1000; % calculate Hertz values

subplot(2,1,1),

plot(f,m); % plot magnitudes

ylabel('Abs. Magnitude'), grid on;

subplot(2,1,2),

plot(f,p*180/pi); % plot phase in degrees

ylabel('Phase [Degrees]'), grid on;

xlabel('Frequency [Hertz]');

The fft()function

calculates the discrete fourier transform in the most efficient way

possible, but transforms of length

2m (m=integer) are considerably fastest

to calculate. Use sizes of 512, 1024, etc for the

fastest speed.

Typically we want the magnitude or the power of a spectrum in

decibels, which can be obtained with:

Y=fft(x,1024);

magnitude = abs(Y);

powerdB=20*log10(magnitude+eps);

To convert back from harmonic amplitudes to a waveform use the

inverse discrete Fourier transform

function ifft(). Be aware that this function takes complex harmonic amplitudes and

deliverscomplex sample

amplitudes. Use the function real() to get the

real part of the signal amplitudes:

Y=fft(x1,1024); % waveform -> spectrum

%%% process Y in some way %%%

x2=real(ifft(Y)); % spectrum -> waveform

3. Windowing a signal

Often we want to analyse a long signal in overlapping short

sections called “windows”. For example we may

want to calculate an average spectrum, or to calculate a

spectrogram. Unfortunately we cannot simply chop

the signal into short pieces because this will cause sharp

discontinuities at the edges of each section. Instead it is preferable to have smooth joins between

sections. Raised cosine windows are a popular

shape for the joins:

You can use the MATLAB function hamming() to

design smooth windows of a given length, and then you can use code

such as the following to divide the signal into sections:

% divide up a signal into windows

x =

sin(2*pi*500*(1:10000)/10000); % example signal

nx =

length(x); % size of signal

w =

hamming(32)'; % hamming window

nw =

length(w); % size of window

pos=1;

while (pos+nw <=

nx) % while enough signal left

y =

x(pos:pos+nw-1).*w; % make window y

%%%% process window y %%%%

pos = pos +

nw/2; % next window

end

4. Spectrograms

MATLAB has a built-in function specgram() for

spectrogram calculation. This function divides a

long signal into windows and performs a fourier transform on each

window, storing complex amplitudes in a table in which the columns

represent time and the rows represent frequency.

Call the specgram function on a signal x with a call:

[B,f,t] = specgram(x,nfft,fs,window,overlap)

Here, nfft is

size of the fourier transform, use this to control the number of

frequency samples on the vertical axis of the

spectrogram; fs is

the sampling rate of the input

signal; window is

the number of input samples per vertical slice, use this to control

the analysis bandwidth; overlap is

the number of samples by which adjacent windows overlap, use this

to control the number of vertical slices per

second. The output B is

a table of complex amplitudes; f is

a vector of frequency values used to label the rows,

and t is

a vector of time values used to label the columns.

Here is a complete example, in which we also use

the imagesc()and colormap() functions

to display the results on a conventional grey scale representing

the top 50dB of the signal:

% get a portion of signal

[y,fs]=wavread('six.wav',[2000 12000]);

%

% calculate the table of amplitudes

[B,f,t]=specgram(y,1024,fs,256,192);

%

% calculate amplitude 50dB down from maximum

bmin=max(max(abs(B)))/300;

%

% plot top 50dB as image

imagesc(t,f,20*log10(max(abs(B),bmin)/bmin));

%

% label plot

axis xy;

xlabel('Time (s)');

ylabel('Frequency (Hz)');

%

% build and use a grey scale

lgrays=zeros(100,3);

for i=1:100

lgrays(i,:)

= 1-i/100;

end

colormap(lgrays);

Reading

An essential introduction to spectral analysis can be found in

Rosen & Howell, Signals and Systems for Speech and

Hearing, Academic Press, 1990; ISBN: 0125972318.

Units 5 and 6 from the Introduction to Digital Signal Processing

course at http://www.phon.ucl.ac.uk/courses/spsci/dsp/ provides

more mathematical detail about fourier analysis and windowing.

The best introduction to signal processing is the book:

Introductory Digital Signal Processing with Computer Applications,

by P.A. Lynn, W. Fuerst, B. Thomas, Paperback - 494 pages 2nd

edition (1998) John Wiley and Sons; ISBN: 0471976318

Exercises

For these exercises, use the editor window to enter your code, and

save your answers to files under your account on the central

server. When you save the files, give them the file extension of

".m". Run your programs from the command

window.

1.Write

a program (ex91.m)

that asks the user for a filename (WAV file), then loads the first

1024 samples, applies a Hamming window and calculates the power

spectrum. Plot the windowed signal and the

spectrum on a single figure with waveform at the top and the

spectrum at the bottom. Make sure you express the

spectral amplitudes in decibels.

2.Adapt

the last exercise (ex92.m)

to calculate and display the average spectrum for the whole

file.

3.Write

a program (ex93.m)

that calculates a high-resolution linear filterbank analysis of a

given file and displays the results on a grey scale, following

spectrogram conventions. Use at least 128

overlapping band-pass filters with bandwidths of 300Hz and an

output rate of 1000 frames/sec.

4.(Homework)

Write a program that filters a signal in the

frequency domain. That is: it takes a signal,

divides it up into overlapping window segments, then for each

segment calculates a FFT, modifies the spectral magnitudes

according to the required frequency response, then performs an

inverse FFT to make a filtered segment. The

filtered segments can then be added together to create the filtered

output signal.

转自：http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect9.html