Randomise phase when using dsp.STFT / ISTFT real time?

Hi,

I would like to create a plugin that performs STFT on the input and then randomises the phases before ISTFT-ing again. I tried to implement this using the dsp.STFT & dsp.ISTFT objects combined with the code for phase randomisation suggested here: https://uk.mathworks.com/matlabcentral/answers/451578-fft-and-ifft-random-phases but it doesn't work (sounds choppy). When I run this I am also getting this message on the command window: Warning: Integer operands are required for colon operator when used as index.

Is there something obvious I am missing and if not is there any other way to implement this ?

I attach a snippet of a real-time processing stream I created to test this, below :

%Create amplitude modulated sine

fs = 44100;
t = 0:1/fs:5;
t = t';
x = (1+cos(2*pi*50*t)).*cos(2*pi*1000*t);
audiowrite('amplmodesine.wav',x,fs)

%---------------------------------------------------
%create real-time audio objects 

FrameLength = 512;
afr = dsp.AudioFileReader('amplmodsine.wav',...
    'SamplesPerFrame',FrameLength);
adw = audioDeviceWriter('SampleRate',afr.SampleRate);

%----------------------------------------------------
%create stft and istft objects

WindowLength = FrameLength;
HopLength = 16;
numHopsPerFrame = FrameLength / 16;
FFTLength = 1024;
win = hamming(WindowLength,'periodic');
stf = dsp.STFT(win,WindowLength-HopLength,FFTLength);
istf = dsp.ISTFT(win,WindowLength-HopLength,1,0);

%-----------------------------------------------------
%processing stream

while ~isDone(afr)
    cleanAudio =  afr();
    y = zeros(FrameLength,1); 
    for index = 1:numHopsPerFrame        
        X = stf(cleanAudio((index-1)*HopLength+1:index*HopLength)); 
        [N,L]=size(X);
% Add random phases
     rnd_theta= -pi+pi.*rand(N,L/2-1); 
     X(:,2:L/2)=X(:,2:L/2).*exp(1i*rnd_theta);
     X(:,L/2+2:L)=X(:,L/2+2:L).*exp(-1i*flip(rnd_theta,2));
% Convert back to time-domain
        y((index-1)*HopLength+1:index*HopLength) = istf(X);        
    end    
% Listen 
    adw(y);
end

Answer 1

[jibrahim]

Hi Angeliki,

I think you need to invert the order here:

[N,L] size(X);

With the STFT object, the second dimension is the number of channels, so L is equal to 1 in your code. I think you want L to be equal to 1024 (fft length) and N to be equal to 1.

The dimensions of rnd_theta would have to change accordingly.

[Angeliki Mourgela]

Hi jibrahim,

Thanks so much for this! It was one of the problems indeed and it now works on the real time stream loop. However, I am still not able to make it work inside the plugin configuration. It sounds very choppy like there are small jumps between frames. Any idea what could be causing this only on the plugin class?

[jibrahim]

Hi Angeliki,

Can you attach your plugin code? I can take a look. It might be dropping audio samples. The audio device writer object returns the number of dropped samples as an output, so you can check that number as you run the plugin in a loop to see if that is happening

[Angeliki Mourgela]

Hi jibrahim,

Sure, here is the code of the plugin. I came across many issues using the previous phase randomisation method, so I was adviced to use the following one. It works on a processing loop but again sounds choppy on a plugin class. I do get 4096 samples underun on testbench.

Let me know what you think!

classdef testplug <audioPlugin
methods
function plugin = testplug()
end
function out = process(plugin,in)   
in = in(:,1);
FrameLength = length(in);
WindowLength = FrameLength;
HopLength = 16;
numHopsPerFrame = FrameLength/HopLength;
FFTLength =  WindowLength;  
win = hamming(WindowLength,'periodic');
stf = dsp.STFT(win,WindowLength-HopLength,FFTLength);
istf = dsp.ISTFT(win,WindowLength-HopLength,1,0);


    
out = zeros(size(in));

for index = 1:numHopsPerFrame        
       
 X_L = stf(in((index-1)*HopLength+1:index*HopLength));
 
 [L,N]=size(X_L);

mag = abs(X_L);%original magnitude
phi = angle(X_L);
rnd_theta = -pi/2+pi/2*rand(L,N);
mix_theta = phi-rnd_theta;
% 
% 
XL_rand = mag.*exp(1i*mix_theta);        
% Convert back to time-domain
out((index-1)*HopLength+1:index*HopLength) = istf(XL_rand);        
end   
end
end
end

[Angeliki Mourgela]

Also, I don't know if it helps but the time scope with a 500 Hz Sine as input, shows this when I run the stft and istft without the phase shifts :

And this with the phase shifts but the original magnitude:

[jibrahim]

Hi Angeliki,

In your plugin code, the STFT and ISTFT objects are recreated every time you call process, so you are not holding the STFT state between process calls. The STFT and ISTFT objects should be properties on the plugin. You can create them in the class constructor. Here is an example:

classdef testplug2 <audioPlugin
    
    properties (Access = private)
        stf
        istf
        win
    end
    
    methods
        function plugin = testplug2()
            WindowLength = 512;
            win =  hamming(WindowLength,'periodic');
            HopLength = 16;
            FFTLength =  WindowLength;
            plugin.stf = dsp.STFT(win,WindowLength-HopLength,FFTLength);
            plugin.istf = dsp.ISTFT(win,WindowLength-HopLength,1,0);
        end
        function out = process(plugin,in)
            
            in = in(:,1);
            FrameLength = length(in);
            WindowLength = FrameLength;
            HopLength = 16;
            numHopsPerFrame = FrameLength/HopLength;
            
            out = zeros(size(in));
            for index = 1:numHopsPerFrame
                
                X_L = plugin.stf(in((index-1)*HopLength+1:index*HopLength));
                
                [L,N]=size(X_L);
                mag = abs(X_L);%original magnitude
                phi = angle(X_L);
                rnd_theta = -pi/2+pi/2*rand(L,N);
                mix_theta = phi-rnd_theta;
                %
                %
                XL_rand = mag.*exp(1i*mix_theta);
                % Convert back to time-domain
                out((index-1)*HopLength+1:index*HopLength) = plugin.istf(XL_rand);
            end
        end
    end
end

[Angeliki Mourgela]

Oh my I'm gonna cry, it actually works now! Thank you ever so much

jibrahim you have been extremely helpful! :)

[jibrahim]

You are welcome!