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);