#!usr/bin/env python # # Programmer: Craig Stuart Sapp # Creation Date: Wed Feb 3 10:22:19 PST 2010 # Last Modified: Thu Mar 4 18:16:56 PST 2010 # Filename: craigspl.py # Syntax: Python 2.6 # """ Full-scale reference amplitude calculations for coordinating the amplitudes of the MDCT and FFT spectra. Source code for answering Question 3 on HW 3 is found in the 'if __name__ == "__main__":' section. """ import numpy as np import numpy.fft as npf from craigwindow import * from craigmdct import * ############################## ## ## calculateDFTreferenceAmp -- Calculate the peak bin amplitude ## when a sinewave is aligned exactly on a frequency bin. Using ## a frequency at 1/2 of the Nyquist frequency to avoid interactions ## with the negative frequency component. ## def calculateDFTreferenceAmp(window, peakamp=1.0): """calculateDFTreferenceAmp -- calculate the corresponding peak amplitude in a DFT spectrum for a given peak sinewave amplitude in the time domain. This amplitude is also dependent on the window type which is the second input. Arguments: window -- array of numbers containing the windowing function. peakamp -- peak amplitude of time-domain sinewave to uas as reference. Returns: the reference amplitude of the peak in the spectrum for the given peak amplitude. """ window = np.array(window) peakamp = float(peakamp) N = len(window) freq = int(N/4) signal = peakamp * np.cos(2.0 * np.pi * freq * np.arange(N) / N) spectrum= npf.fft(signal * window) return max(np.abs(spectrum)) ############################## ## ## calculateMDCTreferenceAmp -- Calculate the peak bin amplitude ## when a sinewave is aligned exactly on a frequency bin. Using ## a frequency as about 1/2 of the Nyquist frequency to avoid ## interactions with the negative frequency component. The MDCT frequency ## bins are shifted 1/2 bin higher than for the DFT, so the frequency ## is on a 1/2 bin in the DFT due to the 1/2 bin shift in the ## freq variable below. ## def calculateMDCTreferenceAmp(window, peakamp=1.0): """calculateMDCTreferenceAmp -- calculate the corresponding peak amplitude in an MDCT spectrum for a given peak sinewave amplitude in the time domain. This amplitude is also dependent on the window type which is the second input. Arguments: window -- array of numbers containing the windowing function. peakamp -- peak amplitude of time-domain sinewave to use as reference. Returns: the reference amplitude of the peak in the spectrum for the given peak amplitude. """ window = np.array(window) peakamp = float(peakamp) N = len(window) freq = int(N/4)-0.5 signal = peakamp * np.cos(2.0 * np.pi * freq * np.arange(N) / N) spectrum= MDCT(signal * window) return max(np.abs(spectrum)) ########################################################################### # # Testing code: # if __name__ == "__main__": """ Homework 3, problem 3 """ import pylab as pl import numpy as np import numpy.fft as npf from craigwindow import * from craigmdct import * N = 1024 freq = 1000.0 srate = 44100.0 peakamp = 1.0 amp = 0.5 sploffset = 96 signal = amp * np.cos(2.0*np.pi*freq*np.arange(N)/srate) recwin = np.ones(N) sinwin = SineWindow(N) hanwin = HanningWindow(N) kb1win = KBDWindow(N, 1.0) kb4win = KBDWindow(N, 4.0) recDFTref = calculateDFTreferenceAmp(recwin, peakamp) sinDFTref = calculateDFTreferenceAmp(sinwin, peakamp) hanDFTref = calculateDFTreferenceAmp(hanwin, peakamp) kb1DFTref = calculateDFTreferenceAmp(kb1win, peakamp) kb4DFTref = calculateDFTreferenceAmp(kb4win, peakamp) recMDCTref = calculateMDCTreferenceAmp(recwin, peakamp) sinMDCTref = calculateMDCTreferenceAmp(sinwin, peakamp) hanMDCTref = calculateMDCTreferenceAmp(hanwin, peakamp) kb1MDCTref = calculateMDCTreferenceAmp(kb1win, peakamp) kb4MDCTref = calculateMDCTreferenceAmp(kb4win, peakamp) ## DFT amplitudes are unormalized (1/2 factor included ## since not doubling the positive frequency amplitudes). # recDFTref # 512.0 # 1 * N/2 # sinDFTref # 325.94945128396904 # 2/pi * N/2 # hanDFTref # 256.0 # 1/2 * N/2 # kb1DFTref # 330.66742639151016 # 0.64583 * N/2 # kb4DFTref # 303.39089615553968 # 0.59256 * N/2 ## MDCT amplitudes are normalized (factor of 2 included in ## transform to double amplitude of positive frequecies) # recMDCTref # 0.99999882345170776 # 1 # sinMDCTref # 0.63661927301775245 # 2/pi # hanMDCTref # 0.49999941172585388 # 1/2 # kb1MDCTref # 0.64583405728406651 # kb4MDCTref # 0.59255964687733975 # ignore negative frequencies in DFT spectra (negative frequencies # have Hermetian symmetry (real:symmetric, imaginary:antisymmetric) absDFTrec = np.split(np.abs(npf.fft(signal*recwin)), 2)[0] absDFTsin = np.split(np.abs(npf.fft(signal*sinwin)), 2)[0] absDFThan = np.split(np.abs(npf.fft(signal*hanwin)), 2)[0] absDFTkb1 = np.split(np.abs(npf.fft(signal*kb1win)), 2)[0] absDFTkb4 = np.split(np.abs(npf.fft(signal*kb4win)), 2)[0] # MDCT output is only positive frequencies (second half of # full transform is anti-symmetric version of first half) absMDCTrec = np.abs(MDCT(signal*recwin)) absMDCTsin = np.abs(MDCT(signal*sinwin)) absMDCThan = np.abs(MDCT(signal*hanwin)) absMDCTkb1 = np.abs(MDCT(signal*kb1win)) absMDCTkb4 = np.abs(MDCT(signal*kb4win)) # convert DFT spectra to DB_SPL splDFTrec = sploffset + 20 * np.log10(absDFTrec/recDFTref) splDFTsin = sploffset + 20 * np.log10(absDFTsin/sinDFTref) splDFThan = sploffset + 20 * np.log10(absDFThan/hanDFTref) splDFTkb1 = sploffset + 20 * np.log10(absDFTkb1/kb1DFTref) splDFTkb4 = sploffset + 20 * np.log10(absDFTkb4/kb4DFTref) # convert MDCT spectra to DB_SPL splMDCTrec = sploffset + 20 * np.log10(absMDCTrec/recMDCTref) splMDCTsin = sploffset + 20 * np.log10(absMDCTsin/sinMDCTref) splMDCThan = sploffset + 20 * np.log10(absMDCThan/hanMDCTref) splMDCTkb1 = sploffset + 20 * np.log10(absMDCTkb1/kb1MDCTref) splMDCTkb4 = sploffset + 20 * np.log10(absMDCTkb4/kb4MDCTref) # Examine dB_SPL maxima in both spectra with various windows. # If the input peak amplitude is 1.0, then the dB_SPL max for # the spectra should be 96 dB (the dbsploffset variable). # However, since the 1000 Hz peak is at bin 23.22, the true # maxima is not directly visible in the sampled spectra. So the # maximum observed amplitudes will be slightly lower (up to 2 dB # lower for sine windowing and 1.4 dB for Hanning windowing in # the worst case). ## For amp = 1.0: ## To get exactly 96, you will have to examine the DTFT (continuous ## form of the DFT) more closely for the 1000 Hz tone at bin 23.22... # max(splDFTrec) # 95.310978661932765 # max(splDFTsin) # 95.604489405903408 # max(splDFThan) # 95.728134825841494 # max(splDFTkb1) # 95.581067925020477 # max(splDFTkb4) # 95.701368021646104 # # max(splMDCTrec) # 94.73811335908681 # max(splMDCTsin) # 95.237247650043884 # max(splMDCThan) # 95.438770142575223 # max(splMDCTkb1) # 95.199062795941018 # max(splMDCTkb4) # 95.394811095218699 ## For amp = 0.5: (6 dB lower than for amp=1.0, so expect 90 dB) ## Notice that an amplitude decrease of 1/2 is equivalent to a decrease ## of 6 dB: 20 * log_10(0.5) = -6.0205999132 dB. ## To get exactly 89.9794, you will have to examine the DTFT (continuous ## form of the DFT) more closely for the 1000 Hz tone at bin 23.22... # max(splDFTrec) # 89.29037874865314 # max(splDFTsin) # 89.583889492623783 # max(splDFThan) # 89.707534912561869 # max(splDFTkb1) # 89.560468011740852 # max(splDFTkb4) # 89.680768108366479 # # max(splMDCTrec) # 88.717513445807185 # max(splMDCTsin) # 89.216647736764259 # max(splMDCThan) # 89.418170229295598 # max(splMDCTkb1) # 89.178462882661393 # max(splMDCTkb4) # 89.374211181939074 # Frequency axis of DFT and MDCT are rotated from each other by 1/2 bin: DFTfreq = srate * np.arange(N/2) / N MDCTfreq = srate * (np.arange(N/2)+0.5) / N ##### ##### ##### Generating plots for HW3, Question 3 ##### ##### ##### ##### Examine Sine/Hanning Windows in DFT/MDCT ##### ##### Plot 3 in Solutions (HW3, Q3) ##### pl.clf() pl.plot( \ DFTfreq, splDFTsin, \ MDCTfreq, splMDCTsin, \ DFTfreq, splDFThan, \ MDCTfreq, splMDCThan, \ ) pl.xlabel("frequency [Hz]") pl.ylabel("Signal Level dB_SPL") pl.title("DFT/MDCT spectra for 1000 kHz sinewave at 0.5 peak amplitude") maxval = int(max(splDFThan+0.5)) pl.text(2000, maxval, 'max = ' + `maxval`) pl.legend(( \ "DFT$_{sin}$", \ "MDCT$_{sin}$", \ "DFT$_{han}$", \ "MDCT$_{han}$", \ ), loc='upper right') pl.grid(True) pl.show() ##### Looking in terms of DFT bins on x-axis: ##### ##### Plot 4 in Solutions (HW3, Q3) ##### DFTbins = np.arange(N/2) MDCTbins = np.arange(N/2) + 0.5 targetbin = freq * N / srate # 23.219954648526077 pl.clf() pl.plot( \ DFTbins, splDFTsin, 'b-', \ DFTbins, splDFTsin, 'bo', \ MDCTbins, splMDCTsin, 'g-', \ MDCTbins, splMDCTsin, 'go', \ DFTbins, splDFThan, 'r-', \ DFTbins, splDFThan, 'ro', \ MDCTbins, splMDCThan, 'c-', \ MDCTbins, splMDCThan, 'co', \ ) pl.xlabel("frequency [bin]") pl.ylabel("Signal Level dB_SPL") pl.title("DFT/MDCT spectra for 1000 kHz sinewave at 0.5 peak amplitude") pl.text(targetbin, 40, `freq`) pl.axvline(targetbin, 0, 100) pl.legend(( \ "DFT$_{sin}$", "", \ "MDCT$_{sin}$", "", \ "DFT$_{han}$", "", \ "MDCT$_{han}$", "", \ ), loc='upper right') pl.grid(True) pl.show() ##### Examine various windows in DFT ##### ##### Plot 1 in Solutions (HW3, Q3) ##### pl.clf() pl.plot( \ DFTfreq, splDFTrec, \ DFTfreq, splDFTkb1, \ DFTfreq, splDFTsin, \ DFTfreq, splDFTkb4, \ DFTfreq, splDFThan, \ ) pl.xlabel("frequency [Hz]") pl.ylabel("Signal Level dB_SPL") pl.title("DFT spectra for 1000 kHz sinewave at 0.5 peak amplitude") maxval = int(max(splDFThan+0.5)) pl.text(2000, maxval, 'max = ' + `maxval`) pl.legend(( \ "DFT$_{rec}$", \ "DFT$_{kb1}$", \ "DFT$_{sin}$", \ "DFT$_{kb4}$", \ "DFT$_{han}$", \ ), loc='upper right') pl.grid(True) pl.show() ##### Examine various windows in MDCT ##### ##### Plot 2 in Solutions (HW3, Q3) ##### pl.clf() pl.plot( \ MDCTfreq, splMDCTrec, \ MDCTfreq, splMDCTkb1, \ MDCTfreq, splMDCTsin, \ MDCTfreq, splMDCTkb4, \ MDCTfreq, splMDCThan, \ ) pl.xlabel("frequency [Hz]") pl.ylabel("Signal Level dB_SPL") pl.title("MDCT spectra for 1000 kHz sinewave at 0.5 peak amplitude") maxval = int(max(splMDCThan+0.5)) pl.text(2000, maxval, 'max = ' + `maxval`) pl.legend(( \ "MDCT$_{rec}$", \ "MDCT$_{kb1}$", \ "MDCT$_{sin}$", \ "MDCT$_{kb4}$", \ "MDCT$_{han}$", \ ), loc='upper right') pl.grid(True) pl.show()