Pimped up audio phase shift demo for https://psychtoolbox.discourse.group/t/how-to-change-the-phase-of-a-sine-wave-with-psychportaudio-in-realtime/4086/23?u=mariokleiner

BasicSoundPhaseShiftDemoPimped.m

function BasicSoundPhaseShiftDemoPimped(showit, targetChannel)
% BasicSoundPhaseShiftDemoPimped([showit=1][, targetChannel=1])
%
% Demonstrates how one can play back a phase-shifted sine tone, with dynamically
% adjustable phase shift, free of audible artifacts during phase shift change.
%
% This uses PsychPortAudio's realtime mixing and some trigonometric math to
% synthesize a cosine wave of selectable phase in realtime from the weighted sum
% of a pair of a cosine + sine wave with different relative amplitudes.
%
% The sine wave is put into one sound channel. A cosine wave of same frequency
% (but 90 degrees phase-shifted) is put in another sound channel. Both sound
% channels are mixed together in realtime by PsychPortAudio, with controllable
% relative volume (== amplitude == mix-weight) for each of the two channels, and
% the result of this two-channel weighted mix is output to the 2nd audio output
% channel of the actual soundcard, creating a synthesized phase-shifted cosine
% wave of same frequency. Changing the relative channel volumes changes the phase
% shift of the cosine output wave in realtime, without artifacts / discontinuities.
% The proper relative volume levels are computed from the desired phase shift and
% set in the internal subfunction updatePhase() by the magic of math.
%
% A second reference sine wave is output to the 1st audio channel of the soundcard,
% if 'targetChannel' == 1, so one can drive two speakers / transducers, one outputting
% a phase-shifted sine wave, relative to th other speaker / transducer. For testing
% purposes, one can set 'targetChannel' == 2 to mix the reference sine-wave with the
% phase-shifted cosine wave both into output speaker channel 2, to simulate potential
% constructive / destructive interference between reference and phase-shifted sine
% wave. Setting 'showit' == 1 will capture the audio signal send to both soundcard
% channel and visualize it in a Psychtoolbox onscreen window, for illustrative purposes
% and for basic debugging. If you want to run a real auditory experiments, you'd use
% 'targetChannel' == 1 and some external microphones and oscillograph to check for
% proper audio output independently of PsychPortAudio, Psychtoolbox and your computer.
%
% One use case for this are studies as described under this link...
% https://psychtoolbox.discourse.group/t/can-i-change-the-phase-of-a-sine-wave-with-psychportaudio/4086
% ...more specifically experiments like this one about air conduction and bone conduction
% of sound:
%
% https://www.researchgate.net/publication/6535047_Simultaneous_cancellation_of_air_and_bone_conduction_tones_at_two_frequencies_Extension_of_the_famous_experiment_by_von_Bekesy
%
% Have a look at BasicAMAndMixScheduleDemo on how to apply amplitude modulation (AM)
% of signals by attaching AM modulator slave devices to the pafixedsine and pashiftsine
% slave devices used here. A suitable time series of AM envelope values would allow to
% gate the output sine waves.
%
% This demo was sponsored by a priority support request. Thanks!
%
% Usage:
%
% ESCAPE key ends the demo.
%
% Left and right cursor keys allow to shift phase interactively.
%
% Optional parameters:
%
% 'showit' If set to 1, visualize actual output signals, channel 1 in red,
%          channel 2 in green. Defaults to 1 for showing output.
%
% 'targetChannel' To which channel should the fixed phase reference signal be
%                 output? 1 = channel 1 (default). 2 = channel 2 will output into
%                 the same channel as the phase-shifted signal, so both signals
%                 will show constructive or destructive interference, depending
%                 on phase-shift of the 2nd signal which always goes to channel 2.

% History:
% 29-Oct-2021 mk  Written.

  if nargin < 1 || isempty(showit)
    showit = 1;
  end

  if nargin < 2 || isempty(targetChannel)
    targetChannel = 1;
  else
    if ~ismember(targetChannel, [1,2])
      error('targetChannel argument must be 1 or 2.');
    end
  end

  % Setup defaults:
  PsychDefaultSetup(2);

  % Define control keys:
  ESCAPE = KbName('ESCAPE');
  leftArrow = KbName('LeftArrow');
  rightArrow = KbName('RightArrow');

  % Perform basic initialization of the sound driver:
  InitializePsychSound(1);

  % Number of physical channels to use on the real soundcard:
  nrchannels = 2;

  % Frequency of waves:
  freq = 500;

% Initial phase shift in degrees of cosine wave:
phase = 0;

  % Open and start physical sound card for playback (1) in master mode (8) = (1+8),
  % with low-latency and high timing precision (1), with auto-selected default []
  % samplingRate for device, to consume and output mixed sound from slaves to
  % 'nrchannels' physical output channels:
  pamaster = PsychPortAudio('Open', [], 1 + 8, 1, [], nrchannels);

  % Retrieve auto-selected samplingRate:
  status = PsychPortAudio('GetStatus', pamaster);
  samplingRate = status.SampleRate;

  % Compute minimum length 'wavedur' seconds of a sound buffer with one or
  % more repetitions of the freq Hz sine wave. If you wanted sound signals
  % of defined length, e.g., also to allow things like applying an AM
  % envelope function to it, by use of AM modulator slave devices, you
  % could just set wavedur to the desired sound duration of the total sound
  % vector. This here is just for memory efficiency...
  %
  % For the minimum duration 'wavedur', make sure to increase wavedur to
  % generate multiple period repetitions of the (co)sine wave ih order to
  % make it fit an integral number of samples, in case one period would
  % need a non-integral number of samples. If nothing else, it may help
  % visualization or debugging / reasoning about it: The if statement is
  % executed, e.g., for combinatios of 'freq' 500 Hz and samplingRate of
  % 44100 samples/sec, where one period of a sine wave would require 88.2
  % samples, ie a non-integral number. Repeating the wave 5x by increasing
  % wavedur * 5, ends up with 88.2 * 5 = 441 samples creating one sound
  % playback buffer with a even number of samples.
  wavedur = 1 / freq;
  
wavedur = 10
  
  nsamples = wavedur * samplingRate;
%  if rem(nsamples, 1)
%    wavedur = wavedur * 1 / rem(nsamples, 1);
%  end

  % Define input vector 'support' for the sin() and cos() functions. Playback of
  % a 'freq' Hz pure sine tone at a sampling rate of 'samplingRate'. We create a
  % waveform of 'wavedur' duration, so that full periods of the sine / cosine fit
  % nicely for inifinitely looped playback, without discontinuities at the beginning
  % and end of the vector. If you also wanted to amplitude modulate the signals with
  % some envelope function, you'd have to create a longer 'support' vector, which
  % repeats the waves for not just one period, but often enough to cover the whole
  % signal duration for the amplitude envelope:
  support = 2 * pi * freq * (0:round(wavedur * samplingRate-1)) / samplingRate;

  if showit
    paoutputcapture = PsychPortAudio('OpenSlave', pamaster, 2 + 64);
    PsychPortAudio('GetAudioData', paoutputcapture, 1);
    PsychPortAudio('Start', paoutputcapture);
    win = PsychImaging('OpenWindow', 0, 0, [0, 0, Screen('Windowsize', 0), 200], [], [], [], [], [], kPsychGUIWindow + kPsychGUIWindowWMPositioned);
  end

  % Create slave device for infinite duration fixed sine tone of freq Hz playback (1).
  % It will output one audio channel (1) via audio channel 'targetChannel' (targetChannel)
  % of the real soundcard:
  pafixedsine = PsychPortAudio('OpenSlave', pamaster, 1, 1, targetChannel);

  % Sine tone, freq Hz:
PsychPortAudio('FillBuffer', pafixedsine, cos(support));

  % Create slave device for infinite duration phase-shiftable cosine tone of freq Hz playback (1).
  % It will output via audio channel 2 of the real soundcard. For this, we let the pashiftsine
  % slave device have two (2) sound channels, both feeding their sound as a mix into the 2nd channel
  % of the pamaster ([2, 2]):
  pashiftsine = PsychPortAudio('OpenSlave', pamaster, 1, 2, [2, 2]);

  % First pashiftsine channel for the mix has a cos() wave, second channel has a sin() wave,
  % 90 degrees phase-shifted. Both waves will be mixed together and output via the 2nd
  % master channel, but with different amplitude / volume. The weighted sum of both 90 degrees
  % shifted waves yields a (co)sine wave of 'freq' Hz, with a phase shifted according to the
  % relative amplitude of both waves. See subfunction updatePhase() for the math behind this:
  srcmixtones = [cos(support) ; sin(support)];
  PsychPortAudio('FillBuffer', pashiftsine, srcmixtones);

  % updatePhase() computes proper relative volumes / amplitudes for the two waves and
  % assign it as channel-volumes to the pashiftsine virtual audio device, so the mix
  % results in a 'phase' shifted cosine wave of 'freq' Hz:
  updatePhase(phase, pashiftsine);

pafixedmodulator = PsychPortAudio('OpenSlave', pafixedsine, 32);
pashiftmodulator = PsychPortAudio('OpenSlave', pashiftsine, 32);

% 'envelope1' sound is wavedur seconds at 0.5 Hz frequency, ie 2 seconds period
% built for a playback device with 'samplingRate' Hz sampling
% rate:
envelope1 = (1 + MakeBeep(0.5, wavedur, samplingRate)) / 2;

% Fill it into standard sound buffer of pafixedmodulator and pashiftmodulator
PsychPortAudio('FillBuffer', pafixedmodulator, envelope1);
PsychPortAudio('FillBuffer', pashiftmodulator, [envelope1; envelope1]);

% Volume control for fixed and shifted sine tone, allows relative adjustment:
PsychPortAudio('Volume', pashiftsine, 0.5);
PsychPortAudio('Volume', pafixedsine, 0.5);

% Master software volume control:
PsychPortAudio('Volume', pamaster, 1);

% "Start" playback of all slaves. This does not actually start them, until the
% pamaster is started, but then they all start in perfect sync:
repeats = 1; % 1 repetition or 0 = infinite
PsychPortAudio('Start', pafixedsine, repeats);
PsychPortAudio('Start', pashiftsine, repeats);
PsychPortAudio('Start', pafixedmodulator, repeats);
PsychPortAudio('Start', pashiftmodulator, repeats);

% Start master, wait (1) for start, return sound onset time in startTime.
% Slaves can be independently controlled in their timing, volume, content thereafter.
% Slaves that are already "started" (but actually are waiting for the pamastr to start),
% will then all start synchronously at startTime:
startTime = PsychPortAudio('Start', pamaster, [], [], 1);

  % Loop for keyboard checking and phase adjustment:
  while 1
    [down, ~, keys] = KbCheck(-1);
    if down
      if keys(ESCAPE)
        break;
      end

      if keys(rightArrow)
        % Phase-shift increase by 1 degree:
        phase = mod(phase + 1, 360);
        updatePhase(phase, pashiftsine);
      end

      if keys(leftArrow)
        % Phase-shift decrease by 1 degree:
        phase = mod(phase - 1, 360);
        updatePhase(phase, pashiftsine);
      end
    end

    if showit
      % Get exactly 0.1 seconds of captured sound:
      recorded = PsychPortAudio('GetAudioData', paoutputcapture, [], 0.1, 0.1);

      % Plot both audio tracks into the window:
      xpos = 1:size(recorded, 2);
      xpos = [xpos xpos(end:-1:1)];
      recorded = [recorded(:, 1:end) recorded(:, end:-1:1)];
      Screen('FramePoly', win, [1 0 0], [xpos' , [100 + recorded(1,:) * 100]']);
      Screen('FramePoly', win, [0 1 0], [xpos' , [100 + recorded(2,:) * 100]']);

      % Show it: Do not sync us to stimulus onset, so we won't slow down to display
      % refresh rate.
      Screen('Flip', win, [], [], 1);
    else
      % Nap a bit, so phase doesn't change that quickly - 50 msecs should do:
      WaitSecs('YieldSecs', 0.050);
    end
  end

  % Stop and close down everything audio related:
  PsychPortAudio('Close');
  if showit
    % Close window:
    sca;
  end

  % Optional plotting code:
  % close all;
  % plot(1:length(recorded), recorded(1,:), 'r', 1:length(recorded), recorded(2,:), 'b');

  % Done, bye bye.
end

function updatePhase(phase, pahandle)
  % Convert phase in degrees into radians:
  shift = phase * pi / 180;

  % Compute relative amplitudes / volumes for both waves, for a resulting
  % wave of amplitude 1.0 and 'phase' degrees (aka shift radians) shift:
  %
  % From https://cpb-us-e1.wpmucdn.com/cobblearning.net/dist/d/1007/files/2013/03/Linear-Combinations-and-Sum-Difference-1xsp3e8.pdf
  % we learn the following trick:
  %
  % Given a1 and a2, there is the following equivalence for the weighted sum of
  % a sine and cosine wave of identical frequency x = 2*pi*freq*t:
  %
  % a1 * cos(x) + a2 * sin(x) = A * cos(x - shift) with
  % A = sqrt(a1^2 + a2^2); and shift = arctan(a2 / a1)
  %
  % so for a desired 'shift' it follows that
  % a2 / a1 = tan(shift)  and as we know that tan(shift) = sin(shift) / cos(shift)
  %
  % therefore:
  %
  % a2 = sin(shift), and a1 = cos(shift), and
  % A = sqrt(sin(shift)^2 + cos(shift)^2) = sqrt(1) = 1, ie. A = 1.
  %
  a1 = cos(shift);
  a2 = sin(shift);

  % Assign new volumes / amplitudes atomically. As channel 1 of pahandle contains
  % a cos(x) wave and channel 2 contains a sin(x) wave, the mixing in PsychPortAudio
  % will create the weighted sum a1 * cos(x) + a2 * sin(x), which is equivalent
  % to an audio signal of cos(x - shift), ie. a cosine wave with 'shift' phase shift:
  PsychPortAudio('Volume', pahandle, [], [a1, a2]);
  fprintf('New phase: %f degrees.\n', phase);
end