working with two variables

3 views (last 30 days)
Tom
Tom on 23 Nov 2011
working with two variables i'm not sure if this is possible, but is it possible to make a 3D plot using an equation based on two variables?
here's my script: -
r_dir=5 %Distance of direct sound path (m)
r_a=0.15 %Radius of rotating drum (m)
d_dir=-5 %Position of centre of rotating drum (on axis) (m)
l_dir=0 %Position of centre of rotating drum off axis/90 degrees from the axis - in the same plane as the rotating drum width ways
t=[0.05] %Time (s)
r_d=5 %Direct distance between receiver and centre of rotating drum (m)
f_r=5 %Frequency/Rotations per second of rotating drum (Hz)
w_r=2*pi*f_r %Angular frequency of the rotating drum
l_f=r_a*sin(w_r*t) %Position of opening on the rotating drum - on the l axis (90 degrees to the axis, so 'off' axis) (m)
d_f=r_a*cos(w_r*t) %Position of opening on the rotating drum - on the d axis (on axis) (m)
r_pd=sqrt((l_f-l_dir).^2+(d_f-d_dir).^2) %Path difference between direct sound path distance and the diffracted sound path distance (m)
r_f=r_a+r_pd %Diffracted sound path distance (m)
c=350 %Speed of sound at 28 degrees celcius with 90% relative humidity
t_y=(r_f-r_dir)/c %Time delay between the arrival of the direct sound and the diffracted sound (sec)
f=20:5:20000 %Frequency of sound from 20 Hz to 20 kHz
w=2*pi*f; %Angular frequency of the sound
AR=r_dir./r_f %Amplitued Ratio - Pressure amplitude of the direct sound divided by the pressure amplitude of the diffracted sound (same as diffracted sound path distance divided by direct sound path distance in metres)
MF=sqrt(AR.^2+1+AR*2*cos(t_y*w))
plot(f,MF)
I picked t=0.05 sec as a constant in the above example, but I'd like to have it so that is a variable too - hence leading to some kind of graph with three dimensions. It would then be a 3D graph (surface) - one axis for the result, one for frequency of sound and one for time.

Accepted Answer

Sven
Sven on 23 Nov 2011
It's not the most efficient code, but I've kept things similar to your original input for clarity. Think now of "MF" as a grid or results, with one pixel for each combination of time t, and frequency f:
f = 20:5:20000; %Frequency of sound from 20 Hz to 20 kHz
w=2*pi*f; %Angular frequency of the sound
t = 0:0.001:0.1; % Set of t locations
MF = zeros(length(t),length(f));
r_dir=5; %Distance of direct sound path (m)
r_a=0.15; %Radius of rotating drum (m)
d_dir=-5; %Position of centre of rotating drum (on axis) (m)
l_dir=0; %Position of centre of rotating drum off axis/90 degrees from the axis - in the same plane as the rotating drum width ways
c=350; %Speed of sound at 28 degrees celcius with 90% relative humidity
r_d=5; %Direct distance between receiver and centre of rotating drum (m)
f_r=5; %Frequency/Rotations per second of rotating drum (Hz)
w_r=2*pi*f_r; %Angular frequency of the rotating drum
for i = 1:length(t)
l_f=r_a*sin(w_r*t(i)); %Position of opening on the rotating drum - on the l axis (90 degrees to the axis, so 'off' axis) (m)
d_f=r_a*cos(w_r*t(i)); %Position of opening on the rotating drum - on the d axis (on axis) (m)
r_pd=sqrt((l_f-l_dir).^2+(d_f-d_dir).^2); %Path difference between direct sound path distance and the diffracted sound path distance (m)
r_f=r_a+r_pd; %Diffracted sound path distance (m)
t_y=(r_f-r_dir)/c; %Time delay between the arrival of the direct sound and the diffracted sound (sec)
AR=r_dir./r_f; %Amplitued Ratio - Pressure amplitude of the direct sound divided by the pressure amplitude of the diffracted sound (same as diffracted sound path distance divided by direct sound path distance in metres)
MF(i,:) = sqrt(AR.^2+1+AR*2*cos(t_y*w));
end
figure
subplot(2,1,1), imagesc(f,t,MF), colorbar, xlabel f, ylabel t
subplot(2,1,2), surf(f,t,MF,'EdgeColor','none'), xlabel f, ylabel t, zlabel result
  2 Comments
Tom
Tom on 23 Nov 2011
Thanks a lot for this. I would be really grateful if you would be willing to help me a little further by just helping me to understand a couple of the coding steps that are new to me. These are: -
row 4 - MF = zeros(length(t),length(f));
row 13 & 21 - 'for' and 'end'
row 13 & 20 - i = 1:length(t) & MF(i,:) = ...
and then the different plot types (I just watched a tutorial video on subplots so I get that bit - it's just the bits after the subplot bit).
If you know some good tutorial videos on these concepts - or just some good sticky posts - that would be perfect.
Many thanks
Tom
Sven
Sven on 23 Nov 2011
Sure, think of it this way... in your original code, your output (MF) was a 1-by-M array of results, where N was the number of frequencies (in f) that you wanted to sample at. In that original code, time (t) was just a single scalar number (ie, 1-by-1). What you then said you wanted to do was repeat your output for *many* different measures of "t". So I also defined "t" as a vector. If you think about how your final output (MF) would look *now*, you'll see that it will be like repeating your original 1-by-M array of results N times (where N is the number of time samples in "t").
The best way to *store* all those results would be in a N-by-M matrix. Each row gives all the frequency results for a given single value of time.
My line 4 command: "MF = zeros(length(t),length(f));"
It simply initialises an empty matrix full of zeros to be the exact size (M-by-N) that we want our final output to be. This is a good MATLAB coding practise when your results will be a big matrix. It's called "preallocation".
Lines 13 and 21. They define a "for loop". Think of looping over every single value for "t" and basically doing your original script (which calculated your MF output for "1" value of t)
Line 13: The "colon" operator. Type "1:20" into the command line, you'll see it makes a vector of all the numbers from 1 to 20.
I think you'll find that the keywords I covered ("for loop", "colon operator", "preallocation" (which is a little more advanced)) are covered well in the "Getting Started" section of the doc. It sounds like you're on the right track though. Hope that's helpful.

Sign in to comment.

More Answers (0)

Categories

Find more on Line Plots in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!