Changing gain in frequency domain creates funny artefacts

Hi,

I need some help with a simulation I am working on in matlab.

I am working on a noise suppression algorithm for speech signals sampled at 
8kHz.

When I modify the amplitude of a number of frequency bins, I get distortions 
in the modified signal when I transfer it back to the time domain.

Here is a copy of my script ....I am hoping that somebody in this group 
could tell me what I am doing wrong.

clc
close all
clear all

%% Load signal
% x: noisy speech signal
% s: clean speech signal
% n: noise signal

addpath '..\audiofiles'

% Load speech signal
TargetFs = 8000;
Tsec = 7;
[s,fs]=wavread('cleanspeech.wav',1);
[s,fs]=wavread('cleanspeech.wav',fs*Tsec);
s = s';
s = resample(s,TargetFs,fs);

% Add white noise to the speech signal
dBSNR = 30;
x = awgn(s,dBSNR,'measured');
n = x-s;


%% Init Variables
% sout: final output signal
sout = [];
% number of samples in a sample block
blockSize = TargetFs/100;
% reset block counter
blockCounter = 0;
% calculate total number of blocks to process
numberOfBlocks = floor(length(x)/blockSize);
% calculate fft size
fftSize = 2^nextpow2(3*blockSize);
% define frame size
frameSize = fftSize;
% init sampleframe for output samples
sampleFrameOut = zeros(1,frameSize);
% init phi angle
phi = 0;

%% Main DSP Loop
while (1)
    audioBlockIn = 
x(1+blockCounter*blockSize:blockSize+blockCounter*blockSize);

    % Do fft of audio block
    audioVector = zeros(1,fftSize);
    audioVector(1:blockSize) = audioBlockIn;
    X = fft(audioVector);

    % simulate gain changes
    gainVec = [abs(cos(phi)) abs(cos(phi)) abs(cos(phi)) abs(cos(phi)) 
abs(cos(phi)) abs(cos(phi)) abs(cos(phi)) ones(1,120)];
    phi = phi + 0.03;

    % Create gain vector
    G = [1 gainVec 1 gainVec(end:-1:1)];

    % Modify spectrum
    Y = G.*X;

    % Overlap-Add
    sampleFrameOut = sampleFrameOut + real(ifft(Y));
    audioBlockOut = sampleFrameOut(1:blockSize);
    sampleFrameOut = [sampleFrameOut(blockSize+1:end) zeros(1,blockSize)];

    sout = [sout audioBlockOut];

    blockCounter = blockCounter + 1;
    if (blockCounter == numberOfBlocks)
        break;
    end

end

wavplay(sout,TargetFs)

Here I do it in the time domain and it sounds fine:

%% Main DSP Loop (TIME DOMAIN)
numberOfFilterCoefficients = 128;
delayLine = zeros(1,numberOfFilterCoefficients-1);
while (1)
    audioBlockIn = 
x(1+blockCounter*blockSize:blockSize+blockCounter*blockSize);

    % Create filter
    h = fir1(numberOfFilterCoefficients-1,0.5+0.4*cos(phi));
    phi = phi + 0.03;

    % Filter audio block
    [audioBlockOut,delayLine] = filter(h,1,audioBlockIn,delayLine);

    sout = [sout audioBlockOut];

    blockCounter = blockCounter + 1;
    if (blockCounter == numberOfBlocks)
        break;
    end

end

Hi Tim,

It works fine when the gain vector does not vary with time. Then there are 
no
artefacts....but as soon as I make the gains time varying (which they would 
have to be
in a noise suppression scheme) I get noise artefacts when I transfer the
modified signal back to the time domain.

On 05/26/2011 11:44 AM, John wrote:

(too bad you took out your context)

> Here I do it in the time domain and it sounds fine:
>
> %% Main DSP Loop (TIME DOMAIN)
> numberOfFilterCoefficients = 128;
> delayLine = zeros(1,numberOfFilterCoefficients-1);
> while (1)
> audioBlockIn =
> x(1+blockCounter*blockSize:blockSize+blockCounter*blockSize);
>
> % Create filter
> h = fir1(numberOfFilterCoefficients-1,0.5+0.4*cos(phi));
> phi = phi + 0.03;
>
> % Filter audio block
> [audioBlockOut,delayLine] = filter(h,1,audioBlockIn,delayLine);
>
> sout = [sout audioBlockOut];
>
> blockCounter = blockCounter + 1;
> if (blockCounter == numberOfBlocks)
> break;
> end
>
> end

I didn't plow through your code, but in general if you multiply in the 
frequency domain that's the same as filtering in the time domain. 
Either your time-domain filter isn't derived from your frequency domain 
filter, or you're suffering from wrapping when you do the ifft.

Try taking the IFFT of _just_ your frequency-domain filter and plotting 
it.  Then take the FFT of _just_ the impulse response of your 
time-domain filter.  Compare the time-domain version of your original 
frequency domain filter, and the frequency-domain version of your 
original time-domain filter.  Ideally, all should match.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" was written for you.
See details at http://www.wescottdesign.com/actfes/actfes.html

Figured it out...

    % Create filter
    h = fir1(fftSize/2-1,0.5+0.4*cos(phi));
    G = fft(h,fftSize);
    phi = phi + 0.03;

    % Modify spectrum
    Y = G.*X;

The problem was that the phase response of the filter also has to be 
included in G. So I can't do it like this:

    % Create filter
    h = fir1(fftSize/2-1,0.5+0.4*cos(phi));
    G = abs(fft(h,fftSize));   % Creates artefacts when taking the abs value 
of the impulse resp.
    phi = phi + 0.03;

    % Modify spectrum
    Y = G.*X;

I thought I could get away with modifying the amplitude of the noisy speech 
signal and then doing the IFFT of the modified spectrum, but obviously I 
also have to change the phase.

Is that correct?

On 05/26/2011 12:23 PM, John wrote:
> Figured it out...
>
> % Create filter
> h = fir1(fftSize/2-1,0.5+0.4*cos(phi));
> G = fft(h,fftSize);
> phi = phi + 0.03;
>
> % Modify spectrum
> Y = G.*X;
>
> The problem was that the phase response of the filter also has to be
> included in G. So I can't do it like this:
>
> % Create filter
> h = fir1(fftSize/2-1,0.5+0.4*cos(phi));
> G = abs(fft(h,fftSize)); % Creates artefacts when taking the abs value
> of the impulse resp.
> phi = phi + 0.03;
>
> % Modify spectrum
> Y = G.*X;
>
> I thought I could get away with modifying the amplitude of the noisy
> speech signal and then doing the IFFT of the modified spectrum, but
> obviously I also have to change the phase.
>
> Is that correct?
>
Looks that way to me.  That 'abs' just hurts to look at, in that context.

Certainly, you have to take phase into account or you'll mess everything up.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" was written for you.
See details at http://www.wescottdesign.com/actfes/actfes.html

On May 26, 3:23&#4294967295;pm, "John" <j...@nospam.thanks> wrote:
> Figured it out...

i can't tell that.

if you're using the FFT and iFFT to do filtering, you have to use
either the overlap-add or overlap-save methods.  i can't tell you're
doing that.  the FIR length must be less than the block size and you
may advance only the difference between the block size and the FIR
length (plus 1 sample, if you want).

normally you will have artifacts when changing the FIR characteristic
between blocks, which is simply a problem of processing block-by-
block.  but if the characteristics do not change greatly between
blocks, you might be able to get away with overlap-add.  i would not
use overlap-save with a time-varying transfer function.

r b-j

> i can't tell that.

well...it seems to work


>if you're using the FFT and iFFT to do filtering, you have to use
>either the overlap-add or overlap-save methods.  i can't tell you're
>doing that.  

What is this then?

    sampleFrameOut = sampleFrameOut + real(ifft(Y));
    audioBlockOut = sampleFrameOut(1:blockSize);
    sampleFrameOut = [sampleFrameOut(blockSize+1:end) zeros(1,blockSize)];

You could change your code such that only the magnitude is changed and
the 'noisy phase' is preserved just as in typical noise reduction.

You may have problems with time-aliasing.
You need to smooth your frequency-domain filter in order to truncate
the length of the corresponding time-domain filter. remember that you
are doing a circular convolution and that constrains the length of
your filter.
I'm not sure it is the problem though as I have only glanced quickly
over your code.
Cheers
niarn

John, you can use a oversampled filter bank for varying with time filtering
in the  frequency domain and to avoid aliasing  artifacts.