Contents

Initial stuff

clear,warning off all
s = tf('s');
% I'll be making a 1.1kHz sinusoidal oscillator, and 2 BPFs going from 50Hz
% to 1kHz and

Make prototype low-pass filter

[z p k]=buttap(2);%Get info for a 2nd order butterworth filter LPF
proto=tf(k*poly(z),poly(p));%Make tf

% Get frequency response info
ww=.01:.1:1e3;
[mag phase]=bode(proto,ww);

disp('Low-pass prototype:'),proto %#ok<NOPTS>

% Display frequency response
figure(1);clf;subplot(211),semilogx(ww,20*log10(mag(:,:)))
title('Bode Diagram'),ylabel('Magnitude (dB)')
axis([0 1e2 -30 10])
subplot(212),semilogx(ww,phase(:,:)/180)
ylabel('Phase (pi radians)'),xlabel('Frequency (rps)')
axis([1e-2 1e2 -1.2 .2])

Low-pass prototype:
 
Transfer function:
        1
-----------------
s^2 + 1.414 s + 1
 

Transmission coefficient and reflection coefficient

t_squared=1/(s^4+1); %Transmission coefficient (2nd order butterworth LPF)
rho_squared=1-t_squared;

pees=pole(rho_squared);
rho=s^2/((s-pees(1))*(s-pees(2)));

disp('rho(s):'),rho %#ok<NOPTS>

% Display reflection coefficient stuff
figure(1);clf;subplot(211)
pzmap(rho_squared)
title('Pole-Zero Map for |\rho(j\omega)|^2')
ax=axis;subplot(212),pzmap(rho)
title('Pole-Zero Map for \rho(s)'),axis(ax)

rho(s):
 
Transfer function:
       s^2
-----------------
s^2 + 1.414 s + 1
 

How it should work in theory - the filters

w1=1;

wca=2*pi*300;
bwa=2*pi*1200;

wcb=2*pi*3500;
bwb=2*pi*4800;

[num den]=lp2bp(k*poly(z),poly(p),wca,bwa);
theory_BPa=tf(num,den);

[num den]=lp2bp(k*poly(z),poly(p),wcb,bwb);
theory_BPb=tf(num,den);

disp('Theoretical Band-Pass Filter:'),theory_BPa %#ok<NOPTS>
disp('Theoretical Band-Pass Filter:'),theory_BPb %#ok<NOPTS>

ww=2*pi*1e1:1:2*pi*1e4;
[maga phasea]=bode(theory_BPa,ww); %#ok<NASGU>
[magb phaseb]=bode(theory_BPb,ww); %#ok<NASGU>

% Display frequency response
figure(1);clf;semilogx(ww/(2*pi),20*log10(maga(:,:)),ww/(2*pi),20*log10(magb(:,:)))
title('Bode Diagram'),ylabel('Magnitude (dB)'),xlabel('Frequency (Hz)')
axis([1e1 1e5 -3 2])
legend('First BPF','Second BPF','Location','NorthEast')

Theoretical Band-Pass Filter:
 
Transfer function:
          5.685e007 s^2 - 4.891e-006 s - 0.0002374
-------------------------------------------------------------
s^4 + 1.066e004 s^3 + 6.396e007 s^2 + 3.789e010 s + 1.262e013
 
Theoretical Band-Pass Filter:
 
Transfer function:
              9.096e008 s^2 - 0.0032 s - 49.35
-------------------------------------------------------------
s^4 + 4.265e004 s^3 + 1.877e009 s^2 + 2.063e013 s + 2.339e017
 

How it should work in theory - the original song

[y, Fs, nbits] = wavread('C:\Documents and Settings\bdieseldorff\My Documents\MATLAB\Rufus_Hallelujah.wav',[ .7e6 1.5e6]); %#ok<NASGU>
x=y(:,1);
soundsc(x,Fs)
dur=length(x)/Fs;

How it should work in theory - after the first bandpass filter

tt=0:1/Fs:dur-(1/Fs);
x2=lsim(theory_BPa,x,tt);
soundsc(x2,Fs)

How it should work in theory - after mixing

x3=x2.*cos(2*pi*1800*tt)';
soundsc(x3,Fs)

How it should work in theory - after the second bandpass filter

x4=2*lsim(theory_BPb,x3,tt);
soundsc(x4,Fs)
wavwrite(x4,Fs,'bp2out_matlab')

Foster's etc leading to the first bandpass (lower)

R1a=160;
% z11=(R1*(1+rho))/(1-rho);
z11=R1a*(2*s^2+s*sqrt(2)+1)/(s*sqrt(2)+1); %Fixed Matlab glitch

% za=z11-s*(sqrt(2))*R1;
za=R1a/(s*sqrt(2)+1); %Fixed Matlab glitch

ya=1/za;

disp('z11(s):'),z11 %#ok<NOPTS>
La=sqrt(2)*R1a;
fprintf('Remove pole at infinity with %0.4g H inductor\n',La)
disp('za(s):'),za %#ok<NOPTS>
disp('Invert')
disp('ya(s):'),ya %#ok<NOPTS>
Ca=sqrt(2)/R1a;R2a=R1a;
fprintf('Implement with %0.3d F capacitor and %0.1g Ohm terminating resistor\n',[Ca R2a])

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/prototype.png','BackgroundColor',bg);
image(a)
axis off;

z11(s):
 
Transfer function:
320 s^2 + 226.3 s + 160
-----------------------
      1.414 s + 1
 
Remove pole at infinity with 226.3 H inductor
za(s):
 
Transfer function:
    160
-----------
1.414 s + 1
 
Invert
ya(s):
 
Transfer function:
1.414 s + 1
-----------
    160
 
Implement with 8.839e-003 F capacitor and 2e+002 Ohm terminating resistor

Transform into the first bandpass (lower)

% The inductor becomes a capacitor in series with an inductor
C1a=bwa/(w1*wca^2*La);L1a=w1*La/bwa;

% The capacitor becomes an inductor and a cap in parallel
L2a=bwa/(w1*wca^2*Ca);C2a=w1*Ca/bwa;

fprintf('We now have R1=%0.3d Ohm, R2=%0.3d Ohm, C1=%0.3d F,\nL1=%0.3d H, C2=%0.3d F, L2=%0.3d H\n',[R1a R2a C1a L1a C2a L2a])
% format this line better

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/bandpass.png','BackgroundColor',bg);
image(a)
axis off;

We now have R1=160 Ohm, R2=160 Ohm, C1=9.378e-006 F,
L1=3.001e-002 H, C2=1.172e-006 F, L2=2.401e-001 H

Simulate grounded inductor for the first bandpass (lower)

% Grounded inductor (L2a) L2a=C_L2*R_L2^2
% uuu---v becomes |Ks:1|---/\/\/---v
C_L2a=2.2e-6;R_L2a=sqrt(L2a/C_L2a);

fprintf('Finally, we get R_L2=%.0g Ohm, C_L2=%0.3g F\n',...
    [R_L2a C_L2a])

% Show Ks:1
bg=[1 1 1];figure(2);clf;
a=imread('Figures/Ks1.png','BackgroundColor',bg);
image(a)
axis off;

Finally, we get R_L2=3e+002 Ohm, C_L2=2.2e-006 F

Foster's etc leading to the second bandpass (higher)

So. LP cuts off at 1kHz. I'll have my multiplier at 1.2 kHz. Then my bandpass should go from 1.1kHz to 2.2kHz

R1b=20;
% z11=(R1b*(1+rho))/(1-rho);
z11=R1b*(2*s^2+s*sqrt(2)+1)/(s*sqrt(2)+1); %Fixed Matlab glitch

% za=z11-s*(sqrt(2))*R1b;
za=R1b/(s*sqrt(2)+1); %Fixed Matlab glitch

ya=1/za;

disp('z11(s):'),z11 %#ok<NOPTS>
Lb=sqrt(2)*R1b;
fprintf('Remove pole at infinity with %0.4g H inductor\n',Lb)
disp('za(s):'),za %#ok<NOPTS>
disp('Invert')
disp('ya(s):'),ya %#ok<NOPTS>
Cb=sqrt(2)/R1b;R2b=R1b;
fprintf('Implement with %0.3d F capacitor and %0.1g Ohm terminating resistor\n',[Cb R2b])

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/prototype.png','BackgroundColor',bg);
image(a)
axis off;

z11(s):
 
Transfer function:
40 s^2 + 28.28 s + 20
---------------------
     1.414 s + 1
 
Remove pole at infinity with 28.28 H inductor
za(s):
 
Transfer function:
    20
-----------
1.414 s + 1
 
Invert
ya(s):
 
Transfer function:
1.414 s + 1
-----------
    20
 
Implement with 7.071e-002 F capacitor and 2e+001 Ohm terminating resistor

Transform into the second bandpass (higher)

% The inductor becomes a capacitor in series with an inductor
C1b=bwb/(w1*wcb^2*Lb);L1b=w1*Lb/bwb;

% The capacitor becomes an inductor and a cap in parallel
L2b=bwb/(w1*wcb^2*Cb);C2b=w1*Cb/bwb;

fprintf('We now have R1=%0.3d Ohm, R2=%0.3d Ohm, C1=%0.3d F,\nL1=%0.3d H, C2=%0.3d F, L2=%0.3d H\n',[R1b R2b C1b L1b C2b L2b])
% format this line better

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/bandpass.png','BackgroundColor',bg);
image(a)
axis off;

We now have R1=020 Ohm, R2=020 Ohm, C1=2.205e-006 F,
L1=9.378e-004 H, C2=2.345e-006 F, L2=8.819e-004 H

Simulating in SPICE

% Main circuit
bg=[1 1 1];figure(1);clf;%subplot('position',[0 0 1 569/700]);
a=imread('Figures/freqshifter.png','BackgroundColor',bg);
image(a)
axis off;

% The hypothetical oscillator I would've built
bg=[1 1 1];figure(2);clf;%subplot('position',[0 0 1 567/1054]);
a=imread('Figures/oscillator.png','BackgroundColor',bg);
image(a)
axis off;

% SPICE simulation agreed very closely with MATLAB simulation although there
% is a small amount of clipping after the first bandpass that didn't happen
% in MATLAB

How I wish it worked - the filters

[z p k]=ellipap(8,.1,60);%Get info for a 2nd order butterworth filter LPF
proto=tf(k*poly(z),poly(p));%Make tf

wca=2*pi*450;
bwa=2*pi*1300;

wcb=2*pi*2500;
bwb=2*pi*1500;

[num den]=lp2bp(k*poly(z),poly(p),wca,bwa);
theory_BPa=tf(num,den);

[num den]=lp2bp(k*poly(z),poly(p),wcb,bwb);
theory_BPb=tf(num,den);

disp('Theoretical Band-Pass Filter:'),theory_BPa %#ok<NOPTS>
disp('Theoretical Band-Pass Filter:'),theory_BPb %#ok<NOPTS>

ww=2*pi*1e1:1:2*pi*1e4;
[maga phasea]=bode(theory_BPa,ww);
[magb phaseb]=bode(theory_BPb,ww);

% Display frequency response
figure(1);clf;semilogx(ww/(2*pi),20*log10(maga(:,:)),ww/(2*pi),20*log10(magb(:,:)))
title('Bode Diagram'),ylabel('Magnitude (dB)'),xlabel('Frequency (Hz)')
axis([1e1 1e5 -30 2])
% legend('First BPF','Second BPF','Location','NorthEast')

Theoretical Band-Pass Filter:
 
Transfer function:
 
0.001 s^16 + 1.875e006 s^14 + 7.21e014 s^12 + 9.929e022 s^10 + 4.603e030 s^8 + 6.345e036 s^6 + 2.945e042 s^4 
                                                                                                             
                                                                             + 4.893e047 s^2 + 1.668e052     
                                                                                                             
-------------------------------------------------------------------------------------------------------------
 s^16 + 1.356e004 s^15 + 3.101e008 s^14 + 2.957e012 s^13 + 3.281e016 s^12 + 2.166e020 s^11 + 1.408e024 s^10 
                                                                                                            
         + 5.945e027 s^9 + 2.177e031 s^8 + 4.753e034 s^7 + 8.999e037 s^6 + 1.107e041 s^5 + 1.34e044 s^4     
                                                                                                            
                                                 + 9.655e046 s^3 + 8.095e049 s^2 + 2.83e052 s + 1.668e055   
                                                                                                            
 
Theoretical Band-Pass Filter:
 
Transfer function:
 
0.001 s^16 + 2.044e-015 s^15 + 4.385e006 s^14 + 3.53e-006 s^13 + 6.394e015 s^12 + 2613 s^11 + 4.331e024 s^10 
                                                                                                             
        + 1.074e012 s^9 + 1.493e033 s^8 + 2.651e020 s^7 + 2.637e041 s^6 + 3.925e028 s^5 + 2.37e049 s^4       
                                                                                                             
                                               + 3.228e036 s^3 + 9.894e056 s^2 + 1.138e044 s + 1.374e064     
                                                                                                             
-------------------------------------------------------------------------------------------------------------
 s^16 + 1.565e004 s^15 + 2.302e009 s^14 + 3.041e013 s^13 + 2.224e018 s^12 + 2.44e022 s^11 + 1.175e027 s^10 
                                                                                                           
         + 1.046e031 s^9 + 3.71e035 s^8 + 2.58e039 s^7 + 7.155e043 s^6 + 3.665e047 s^5 + 8.243e051 s^4     
                                                                                                           
                                                + 2.781e055 s^3 + 5.194e059 s^2 + 8.713e062 s + 1.374e067  
                                                                                                           
 

How I wish it worked - the original song

[y, Fs, nbits] = wavread('C:\Documents and Settings\bdieseldorff\My Documents\MATLAB\Rufus_Hallelujah.wav',[ .7e6 1.5e6]);
x=y(:,1);
soundsc(x,Fs)
dur=length(x)/Fs;

Looking at the Frequency content

win=hamming(2^12);
out=stft(x,win,length(win));

freqs=(Fs/length(win))*(1:length(out(:,1)));
figure(2);clf;loglog(freqs,2*abs(out)/length(win))
xlabel('Frequency (Hz)')

How I wish it worked - after the first bandpass filter

tt=0:1/Fs:dur-(1/Fs);
x2=lsim(theory_BPa,x,tt);
soundsc(x2,Fs)

Looking at the Frequency content

win=hamming(2^12);
out=stft(x2,win,length(win));

freqs=(Fs/length(win))*(1:length(out(:,1)));
figure(2);clf;loglog(freqs,2*abs(out)/length(win))
xlabel('Frequency (Hz)')

How I wish it worked - after mixing

x3=x2.*cos(2*pi*1800*tt)';
soundsc(x3,Fs)

Looking at the Frequency content

win=hamming(2^12);
out=stft(x3,win,length(win));

freqs=(Fs/length(win))*(1:length(out(:,1)));
figure(2);clf;loglog(freqs,2*abs(out)/length(win))
xlabel('Frequency (Hz)')

How I wish it worked - after the second bandpass filter

x4=2*lsim(theory_BPb,x3,tt);
soundsc(x4,Fs)
wavwrite(x4,Fs,'bp2out_wishful')

Looking at the Frequency content

win=hamming(2^12);
out=stft(x4,win,length(win));

freqs=(Fs/length(win))*(1:length(out(:,1)));
figure(2);clf;loglog(freqs,2*abs(out)/length(win))
xlabel('Frequency (Hz)')

Conclusion

% I really like the pictures of the frequency content. It seems like the
% sound should be really good. Except that, in the time domain, it really
% isn't. That filter is really sharp. I am, coincidentally, really
% disappointed. Grrr... it seems I was misguided to start with. I think
% losing bandwidth early on makes the sound terrible later. I guess if I
% wanted to actually do this, I'd have to pitch shift it up extremely high
% and then shift it partway back and use awesome fiters everywhere.

Corrected Conclusion

% It actually works well. I'm just not doing what I intended to do. I'm
% doing frequency shifting instead of pitch shifting. But yeah awesome.
% Listen to some wav files.

Clean up

warning on all
clear

Published with MATLAB® 7.7