POWER SPECTRAL DENSITIES AND SPECTROGRAMS OF NOISE BANDS AT FS = 4000 HZ

Contents

 Created by Juan-Pablo Caceres 21-Jul-2006.
 Last Revision: Under Subversion Control System
clear all
close all

Motivation

This document compares the Power Spectral Densities of Noise bands at different frequencies generates using fs = 4000 Hz. It uses ERB to estrimate the critical bandwith. It shows that at this fs, the design has the same maximum power for all the bands. It alse explains some issues concerning the playback of this bands at this small sampling rate.

fs = 4000; %Hz

Center Frequencies

fc1 = 100; fc2 = 300; fc3 = 800;

fm = 0; amp = 0;
dur = 5; %seconds

Center frequency fc1

delf = erb(fc1);
[b,a] = bpdiir(fc1,delf,fs); %design IIR filter design
noisefm1 = noisefm(fs,b,a,dur,fm,amp); %FM modulation
clear b a;

Center frequency fc2

delf = erb(fc2);
[b,a] = bpdiir(fc2,delf,fs); %design IIR filter design
noisefm2 = noisefm(fs,b,a,dur,fm,amp); %FM modulation
clear b a;

Center frequency fc3

delf = erb(fc3);
[b,a] = bpdiir(fc3,delf,fs); %design IIR filter design
%freqz(b,a)
noisefm3 = noisefm(fs,b,a,dur,fm,amp); %FM modulation
clear b a;

Power Spectral Densities Estimates

[Pxx1,w] = pwelch(noisefm1);
[Pxx2,w] = pwelch(noisefm2);
[Pxx3,w] = pwelch(noisefm3);

figure
plot(w/pi*fs/2,10*log10(Pxx1),'b');
grid
xlabel('Frequency (Hz)'); ylabel('Power/frequency (dB/rad/sample)');
title('Power Spectral Density Estimate via Welch');
hold on;
plot(w/pi*fs/2,10*log10(Pxx2),'g');
plot(w/pi*fs/2,10*log10(Pxx3),'r');
legend(['PSD fc = ' num2str(fc1) 'Hz'],...
       ['PSD fc = ' num2str(fc2) 'Hz'],...
       ['PSD fc = ' num2str(fc3) 'Hz']);

Spectrogram

figure
spectrogram(noisefm1+noisefm2+noisefm3,2^9,[],2^12,fs,'yaxis');
xlabel('Time (sec)');
title('Spectrogram of the Noise Bands')

NOTE Concerning Playback at 4000Hz

At least in my computer soundcard, I get aliasing when I playback the noises at 4000Hz. This is an artifact from the soundcard that is not present in the generated noise. To avoid this, before listening, resample the signal using matlab's resample function, e.g.:

noisefm3_44100 = resample(noisefm3,44100,fs,40);