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

Question

Angeliki Mourgela 2020-7-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/565604-randomise-phase-when-using-dsp-stft-istft-real-time

评论： jibrahim 2020-7-23

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

jibrahim 2020-7-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/565604-randomise-phase-when-using-dsp-stft-istft-real-time#answer_468474

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Angeliki Mourgela 2020-7-21

在 MATLAB Online 中打开

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

jibrahim 2020-7-23

在 MATLAB Online 中打开

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 2020-7-23

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

jibrahim you have been extremely helpful! :)

jibrahim 2020-7-23

You are welcome!

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论