How to Generate Band Limited Gaussian White Noise in MATLAB? In some papers, I read that the additive noise is band limited Gaussian white.
How can I simulate this type of noise use MATLAB?

You would generate bandlimited Gaussian noise by first generating white noise, then filtering it to the bandwidth that you desire. As an example:

% design FIR filter to filter noise to half of Nyquist rate
b = fir1(64, 0.5);
% generate Gaussian (normally-distributed) white noise
n = randn(1e4, 1);
% apply to filter to yield bandlimited noise
nb = filter(b,1,n);

Producing full spectrum white noise and then filtering it is like you want to paint a wall of your house white, so you decide to paint the whole house white and then paint back all house except the wall. Is idiotic. (But has sense in electronics).

I made a small C program that can generate white noise at any frequency and any bandwidth (let's say at 16kHz central frequency and 2 kHz "wide"). No filtering involved.

What I did is simple:
inside the main (infinite) loop I generate a sinusoid at center frequency +/- a random number between -half bandwidth and +halfbandwidth, then I keep that frequency for an arbitrary number of samples (granularity) and this is the result:

White noise 2kHz wide at 16kHz center frequency

Pseudo code:

while (true)
{

    f = center frequency
    r = random number between -half of bandwidth and + half of bandwidth

<secondary loop (for managing "granularity")>

    for x = 0 to 8 (or 16 or 32....)

    {

        [generate sine Nth value at frequency f+r]

        output = generated Nth value
    }

}


Finding a signals complimentary sound / frequency I know about complimentary colors and how to get them, since colors are basically frequency / wavelengths.  Is it possible to find a signals complimentary, additive or subtractive frequencies in an audio file / from a signal.  I know how to use FFT to get the phase, amplitude, and frequency of a signal.  I'm just looking for a general idea of how I would go about creating similar color charts but for sound, we have a color charts of complimentary, additive or subtractive colors do we have one for the sound / audio range of humans?

Yes and no. Only monochromatic colors consist of a single wavelength, and those do not really "wrap around" in the form you describe. Color addition and subtraction have more to do with the three types of cone cells in our eyes than with wavelengths. Since there is no audio analogue to these cone cells, the concept of complimentary colors doesn't exist in the same way.

However, if you consider notes instead of frequencies, in the sense that two frequencies that are one or more octaves apart constitute the same note, then you can say that two notes are complimentary if they are a tritone apart.

A couple of ideas:

Bass, mid, treble
Transients, noise, sinusoids


Time period to Instantaneous frequency I have simple smooth time signal with different time period. I would like to transfer period over the time to frequency, and to track frequency changing over the time, I would call it Instantaneous frequency?

PF-period frequency (1/T)... any suggestion?

You can compute the instantaneous frequency of your data in Matlab as follows:

z = hilbert(x);
N = length(z);
f = angle(z(2:N).*conj(z(1:N-1)))/2/pi;

To get units of Hertz, multiply f by the sample rate.

There are other, similar estimators that do a similar thing. This matlab code can help.

Eric Jacobsen also has a couple of pages here and here.

I think you want to create some sort of spectrogram. The only difference is that the spectrogram contains a full frequency spectrum for every time period, and you want to get just the dominant frequency, but you can run a dominant frequency search (maximum search), and you will get what you want. You have control both of the time and frequency resolution.


How to initialize an IIR filter in Matlab I designed an IIR filter and want to do zero, step and projection initialization.  There is mention in different reference papers but no explicit examples for Matlab.  How do you do IIR initialization, for the previously mentioned types, in Matlab?

Not easy. The Matlab filter() function implements a general purpose IIR filter and allows you to pass in an initial state. However, they don't publish the exact definition of the state and the algorithm deployed. It seems to be close to a transposed form II filter but with some subtle numerical differences, so the state doesn't quite match.

In any case it would be better to break this down into second order sections using sosfilter(). However that doesn't allow you to set an initial state. You may have to write your own filter function for that.

YOu still can use filter() to some extent. The default is zero initialization. You can also do step initialization by simply running a constant input for a while and collecting the state as an output argument from the filter, like this:

[b a] = butter(4,1000*2/44100);
% create state for a constant input of "1"
[~,state1] = filter(b,a,ones(1000,1));
% now filter for a constant of 2 but with a state for "1" input
ystep = filter(b,a,2*ones(100,1),state1);
plot(ystep);

Not sure about projection initialization, though.

For your own implementation see:

http://www.emt.tugraz.at/publications/diplomarbeiten/da_hoebenreich/node21.html#fig:StateSpaceDiagram


How Can I Resample a Signal with an Arbitrary Factor (For Example - 128000Hz to 16000.1Hz) in MATLAB? I need to simulate the sampling of a continuous (fsCtu = 128000Hz), acoustic signal with two microphones that have a slight offset in sampling rate (fsMic1 = 16000, fsMic2 = 16000.1) in Matlab. What is the best way to do this?

Things I tried:

The Matlab "resample" command only works for resampling to 16000, not to 16000.1
"interp1" doesn't seem to be an option because I think I need to use bandlimited interpolation for a correct simulation. (Is this assumption correct?)
I tried to write my signal to a wav file and resample it via a system call using this software, then load the processed file with wavread. I'm not sure if this is a good solution. A quick test revealed that this method doesn't give the same result as the "resample" command for resampling to 16000Hz, which I find strange.

Any ideas or suggestions?

There are a couple of ways that you can do it.  The first is with resample, but it is a multi-step process.  First, you have to figure out which interpolation and decimation factors will get you the sample rate you want.

[n, k] = rat(16000.1 / 128000);

That gets you an interpolation factor of 20000 and a decimation factor of 159999.  You factor those to break them up into smaller chunks.

nFactors = factor(n)
kFactors = factor(k)

It turns out that n factors to $2^5 * 5^4$ and k factors to $3 * 7 * 19 * 401$.  All of those are doable, though the decimation by 401 will not have great filtering properties.  Anyway, if you resample in stages you can get the final sample rate you want.

The other way to do it is polynomial interpolation.  Essentially you model your signal as a polynomial through curve fitting techniques, and then you can simply feed in the time values that you want and out will pop the signal values.  This technique can be very effective, but modelling the signal well is a bit of an art.  In particular you don't want to overfit.

The correct way to do this in Matlab is:

128000 -> 16000 ==> resample(x,1,8)

128000 -> 16001 ==> resample(x,16001,128000)

Having said that, I don't know what limitations Matlab places on P & Q.

John


Averaging magnitude squared coherence across multiple time series In a previous post, A. Donda had suggested that, in order to calculate the average magnitude squared coherence of more than one pair of time series (e.g. y1 and x1, and y2 and x2), one ought to follow the following steps:

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);

sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);

sxy2 = cpsd(x2, y2, w, o, n);

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;

sxya = (sxy1 + sxy2) / 2; % NOTE THIS LINE %
ca = abs(sxya) .^ 2 ./ sxa ./ sya;

...where w,o,n are the window parameters, and the other variables should (hopefully) be self-explanatory.

Let me illustrate why this method can be problematic with an example, in the hope that someone might be able to help me with a solution.

Imagine two pairs of signals y1, x1, and y2,x2. All four signals contain the same frequency component (e.g. a 440Hz sine wave). In one pair, the sines waves are in phase w.r.t. each other. In the other pair, the sine waves are in inverse polarity.

Separate MSCs of both of these pairs of signals will yield approx. '1' at 440Hz in both results, identically, telling us that 440Hz is common to both signals in both pairs (i.e. the relative phase of each of the signals within a pair is ignored). If we now wish to calculate the average MSC, we must first compute the average PSDs and then the average CPSDs. Ideally, this average MSC for the two pairs would also yield a '1' at 440Hz. BUT the CPSD output is complex, so averaging the two CPSDs of the two pairs of signals according to the method above will result in cancellation, because the averaging is taken BEFORE the magnitude of the CPSDs is calculated.

I have offered one solution below - is this correct?

I believe I might have a solution... Editing A Donda's work above, I think a simple rearrangement of the last two lines is necessary:

sxya = (abs(sxy1) + abs(sxy2)) / 2;
ca = sxya.^ 2 ./ sxa ./ sya;

So here, the average of the magnitudes is taken, rather than the average of the complex CPSDs. Then, this average can be squared.


OpenCV: Problem at creating a Skeleton of Lines I need to get the lines from a photographie of a sheet of paper. While recherching, I came along this, where one answer referred to this blog entry.
Now I am having trouble with this approach. I have implemented it as follows (with the c++-interface of OpenCV:

// the image has to be grayscale
if (image.channels() != 1)
{
    cv::cvtColor(image, image, cv::COLOR_BGR2GRAY);
}

// we need to enhance the lighting before we can threshold the photography
cv::equalizeHist(image, image);
// a binary image is needed
cv::threshold(image, image, threshold, 255, cv::THRESH_BINARY);
cv::imwrite("thresholdedImage.jpg", image);

// the resulting skeleton
cv::Mat skeleton(image.size(), CV_8UC1, cv::Scalar(0,0,0));
// needed if in-place processing is not possible
cv::Mat temp(image.size(), CV_8UC1);
// eroded image is saved here
cv::Mat eroded;
// needed for morphological transforms (erodation, dilation)
cv::Mat element = cv::getStructuringElement(cv::MORPH_CROSS, cv::Size(3,3));

bool done = true;
int i = 0;
while(done)
{
    // eroding + dilating = opening
    cv::erode(image, eroded, element);
    cv::dilate(eroded, temp, element);
    cv::subtract(image, temp, temp);
    cv::bitwise_or(skeleton, temp, skeleton);
    eroded.copyTo(image);

    done = (cv::countNonZero(image) == 0);
    ++i;
}

For proof of concept, I decide to use a computer-made drawing for testing purposes, to eliminate issues with lighting at this early stage. This can be seen here:

The threshold is correct, but unfortunately the skeleton looks very strange:

I have drawn the skeleton after 100 iterations, where it looked like this:

Do you have any suggestions?

You are actually drawing the skeleton of the background (brighter region).

change

cv::threshold(image, image, threshold, 255, cv::THRESH_BINARY);

in your codes to

cv::threshold(image, image, threshold, 255, cv::THRESH_BINARY_INV);

You should be able to get the right skeleton.

You could well refer here for a few different implementations, which provide acceptable results.


KLT tracker implementation I'm working on a object tracker based on KLT algorithm using matlab.
Can anybody guide me how I can find the features and calculate the distance in the sequence of the frames?
Thanks.

Try this example of tracking a single face using KLT.  Try this for tracking multiple faces.


How do I compare 2 audio files, returning a percentage of difference Is there a way to compare 2 audio files, whereby a percentage of difference will be returned?
I am using C# VS08, .net framework 3.5.

Edited: difference in the sound (like for eg. Audio 1: "HELP", Audio 2:"HELP ME PLEASE", it will return the difference in percentage between this 2 audio.)

I am currently saving recorded files in wav format, and its bit rate and quality will also be the same as what will be recorded and compared against. I wish to ask if there is a way to compare the sound (not just the binaries) in 2 audio files of the same properties and same format (WAV), and show the percentage difference between the 2 sounds in the 2 audio files.

Help will be greatly appreciated.

Could you please define the "percentage of difference"?

For example, please tell us which value do you expect to get:

Between two versions of the same file encoded with different settings?
Between two versions of the same file, except that one is slightly faster than the other?
Between two recordings of the same audio input, but through different microphone/recording equipment?
Between two recordings of the same person saying the same word?
Between two recordings of the same person saying the same word, with a different prosody (rhythm/melody)?
Between two recordings of different persons saying the same word?
Between two recordings of the same person saying the same sentence except for a few words?
Between a recording of a person and that of a cow?
Between a recording of an airplane and music?

What do you want to measure? The speech contents (words)? The melody, rhythm? The overall audio similarity? Encoding/recording equipment differences?

I am using one javascript function for compare tow audio file. By using the same logic you are able to compare any file in any language.

function compireAudio(){
var audio1 = "http://soundjax.com/reddo/86502%5Ealarm.mp3";
var audio2 = "http://soundjax.com/reddo/44368%5EALARME.mp3";
var i,j,d;
var matching = 0;
var t = 0;var i,j,d;
var matching = 0;
var t = 0;
var audio1Arr = Array();
var audio1Len = audio1.length;
for (i = 1; i<=audio1Len; i++)
{
    //reverse so its like a stack
    d = audio1.charCodeAt(audio1Len-i);
    for (j = 0; j < 8; j++)
    {
        audio1Arr.push(d%2);
        d = Math.floor(d/2);
    }
}
var audio2Len = audio2.length;
for (i = 1; i<=audio2Len; i++)
{
    //reverse so its like a stack
    d = audio2.charCodeAt(audio2Len-i);
    for (j = 0; j < 8; j++)
    {
        if(d%2 == audio1Arr[t])
        {
            matching++;
        }
        d = Math.floor(d/2);
        t++;
    }
}
var avarage = Number(matching)/((Number(t)+Number(audio1Arr.length))/Number(2))*Number(100);
alert('The Matching with the two audio is '+avarage+' %.');

}


Applying a high-pass filter on constant signal I recently implemented a Butterworth high-pass filter (2nd order). Everything seems to work fine, except that I have a question when applying such a filter to a constant signal. As there are no high-frequency components in a constant signal, I would expect the filter to yield a constant 0 signal.

In the plots below are my results. There seems to be ``ripple'' in the first couple of frames before the high-pass filtered signal (correctly) converges to 0.

Is this a logical result from applying a Butterworth high-pass filter on a constant signal, or might there be a bug in my code?

EDIT: I've created the same filter in Octave, resulting in the same output. This indicates that we're indeed looking at the step response, as Paul R. indicates.

Octave output

The input to your filter is essentially a step function (since it has value 0 prior to t = 0, and a positive value for t > 0), so you see the step response of the filter, hence the initial ringing. This is expected behaviour, and after a suitable amount of time the step response will have settled to zero.

You used to low-order filter(2nd order). This can be reason that composing high ripples at starting frequency points.I think you can increase the order of filter and try it on your design.


Polyphonic detection/Mulit-pitch detection/Chord recognition I build application witch will be transform piano sound to musical notes. I found some pitch detection algorithms with implementations (mainly based on autocorrelation or fft), but chord recognition, of course is  not working.
What is best method for instrument as piano(both hands + chords)?
Is there any implementation of this problem?

As @jan pointed out, you're probably asking for a little too much, especially if you're looking for a ready implementation. Doing a quick search on Google, I came across several papers that may be a helpful start.

In this paper called Multi-pitch Detection Algorithm Using Constrained Gaussian Mixture, the authors use the Expectation Maximization algorithm to solve a Gaussian Mixture Model to detect multiple tones. This is most likely not a very computationally efficient algorithm, since mixture models are normally very hard to solve, and have to be primarily done offline (no real-time version).

If you're looking for a moderately robust, but potentially fast algorithm, look at this Stack Overflow answer. The author recommends MUSIC and ESPRIT algorithms.

The best thing to do is to give something a try and come back to us with specific DSP and algorithms questions and sites like Stack Overflow for specific implementation questions.

Have you looked at this paper?

https://arxiv.org/pdf/1804.02918&ved=2ahUKEwj46ZyDrajjAhXkIzQIHdt5AqYQFjACegQIAhAB&usg=AOvVaw3OoDyPWKVPW7fDitFpn_cH&cshid=1562693079932


Inter aural time difference generation using sine waves I have a stereo sound for Inter  aural time difference consisting of sine wave(16 samples as 16 bit PCM Data) for left and right audio. From time to time I would like to advance either the left or right audio so that sound in one ear comes earlier than the other. If I just put another sine wave of different frequency in between my normal sine samples to give an illusion of changing sound in one ear, I hear a click. I am playing a buffer of 1024 samples at a time containing looped sine wave samples. If I drop one sample it generates a click as well. Can someone suggest a way of speeding up the Sine wave without generating a click.
e.g loop the sine wave separately for both ears. Now advance the sine wave in one ear so that it feels like a frequency change and then play the normal sound again in both ears. I have fixed sine samples with me. (I work in java so can't do any heavy computations). I have tried to clarify but the question could still be confusing. Ask for specific details.

Thanks

I understand the point.
If you are outputting signal[i] you cant jump to signal[i+n] instantaneously.
You should just fade it to the new point with some window function. Allow some period of fading time.

For example: If I currently output signal[i] and I want to shift it to signal[i+n] in 1 second. I will use half of Hann window to fade it in and out. So I time scaling Hann from pi to 1 second.
and output $signal[i+n]*Hann[i]+signal[i]*(1-Hann[i])$ for 1 second before I switch to signal[i+n] at full scale. This should not generate a click.

You may want to try other window function or different fading period.


What do we call the areas with "1" in a binary image I am looking for possible professional terms. In a binary (1 or 0) image, normally, the areas with "1"s are where the information is. I wonder if there is any professional name for it in image processing. Thank you!

This problem born when you use pixel intensity as physical mass, for example to calculate mass center or inertia moments, in this case the background should be 0 and object should be 1.

Btw your question is not clear because in a binary image {0,1} have the same information content and generally there are not difference.

There is not a name for 1 pixel and 0 pixel, often in common convention 1 is white and 0 is black, but it's not a general rule. Maybe you can call 1 high value pixel and 0 low value pixel but I think that are not common names and they can make confusion in some contexts.

A binary image conventionally distinguishes "foreground" from "background", or "objects" from "background".

More often than not, the foreground pixels, which form the objects, have a value of 1, but some libraries will use a value of 255 (storing the binary image in a 8-bit unsigned integer array). Conversely, the background typically has a value of 0.

The reason to select 0 for background and a non-zero value for foreground is that (at least in every programming language I know) 0 is evaluated as false and non-zero as true. This makes subsequent logic more intuitive.

I have personally never come across a scientific paper on image processing or analysis that reverses this convention, but I have seen plenty of papers that print background as white and foreground as black -- When most pixels are background, this makes the foreground stand out better in print.

Nonetheless, I'm sure there are programs out there that use 0 for foreground. It's not written in stone. It is simply the very common convention.

(Edit to preserve the comment I wrote to the other answer.)

Note that these terms and conventions go back to the beginnings of the field of digital image processing, and probably earlier too.

Azriel Rosenfeld was one of the pioneers of the field. In this book from 1976,
"Digital Picture Processing" (Ch. 8), he says

  [...] the function which has the value 1 for points belonging to the object, and value 0 elsewhere.

And a little later

  [...] deciding whether a given point of the picture belongs to an object or to the background.


ROC curve with classification problems I am working with face anti-spoofing. For classification, the Euclidean distance measure is used. The test image vector is compared with all real image vectors , and the minimum distance is chosen (min.1) . The same vector is compared with all fake image vectors and the minimum value is chosen (min. 2). The minimum rule is applied again to choose the minimum of both min.1 and min.2 to identify the correct class of the test image. My problem is in plotting the ROC curve. I do not have a threshold which can be changed. I am using hard decision rules. How can I plot the curve in this case

Yes. You can.
ROC curve is a Pd vs Pfa plot where Pd and Pfa are the probability of detection and probability of false alarm respectively.


How to average MFCCs from different speakers saying the same word? I am trying to create a simple recognition system. i have recordings of several speakers saying the same word, let's say "west". I compute the MFCC feature vectors of word "west" from each speaker. If i have 10 speakers, i then have 10 matrices, each of them of dimensions N_frames * 13 (13 coeff and N_frames varies from speaker to speaker). I would now like to store a final matrix of dimension
M * 13 which will serve as a reference pattern for word "west". My question is, how to "average" or combine MFCCs from different speakers? (i call this process some kind of training)

What you can try is the DTW (Dynamic Time Warping) distance (https://en.wikipedia.org/wiki/Dynamic_time_warping).

for each instance of the spoken word "west", you compute MFCC's per frame and then calculate the DTW distance between each of the pre-recorded instances of word. this gives the rough estimate of the valid DTW distance. now when a new instance of the word comes we again calculate DTW with one of the earlier instance and if the DTW distance lies roughly in the range calculated before, it is the same word.


Does clipping noise colour the spectrum Does extreme clipping white noise colour the audio spectrum?

I thought (and hoped!) not on the basis of this from research.net

  Clipping is mathematically equivalent to the addition of impulse
  noise, (at the point of clipping you are effectively adding a negative
  impulse to the original waveform); consequently the effect on the
  spectrum is similar to additive impulse noise, with a fairly uniformly
  flat spectrum (depending on the degree of clipping). The spectrum is
  obviously going to be coloured if the clipping occurs with some
  periodicity

Emphasis mine.

Clipping is a non-linear operation and it introduces new frequencies, so it will change the spectrum's shape. The text you quoted is pretty clear on that.

Just as a quick demo, I generated flat noise in Matlab, limited in frequency to the 100 to 300 Hz range. The spectrum is:

Note that the spectrum is close to zero outside the specified range. Also, this noise is not white; it is flat (or would be if I could simulation an infinite number of samples). White noise, by definition, covers the entire frequency range.

I clipped all amplitudes larger than 0.75 to 0.75, and the resulting spectrum is:

I think the difference is clear.

White noise has several definitions.  If you take the definition where all samples are independent, clipping will change the distribution of the individual samples but not make them dependent.  The result will again be independent samples and consequently white noise according to that definition.


Artifact in Bicubic interpolation I wrote an algorithm to do bicubic interpolation of an image.
I used the method desribed in the wikipedia page.
On simple images, the result looks good, but on more complex ones, I got strange artifacts in non-smooth zones.

What could be a probable source of those artifacts, and what can I try to get rid of them ?

My code is on github, but I checked it several times, so I don't think the problem comes from it but rather from the simple implementation suggested in wikipedia

I'm using this particular equation (3rd-order Hermite polynomial interpolation) to do the interpolation (in one dimension):

$$x(n+t)=\frac{1}{2}
\begin{bmatrix}
1 & t & t^2 & t^3 \\
\end{bmatrix}
\begin{bmatrix}
0 & 2 & 0 & 0 \\
-1 & 0 & 1 & 0 \\
2 & -5 & 4 & -1 \\
-1 & 3 & -3 & 1 \\
\end{bmatrix}
\begin{bmatrix}
x[n-1] \\
x[n] \\
x[n+1] \\
x[n+2] \\
\end{bmatrix}$$

where $x[n-1]$, ... $x[n+2]$, .. are the values of the surrounding pixels (in one dimension), and $t$ is the fractional coordinate (where $0 \le t < 1$) of the interpolated pixel between the adjacent pixels $x[n]$, $x[n+1]$. $n$ is an integer and $x(n)=x[n]$ for any integer $n$.

The expected output:

My result, look closely at the bird head,a lot of white pixels have appeared.

---EDIT---
To show that my code does produce coherent results, here is the example of my scaling up on simple image with 16 red pixels.

And here is the ouput of the library (which also use 3rd-order Hermite polynomial interpolation):

And here is my results:

In this case the difference are IMO totally acceptable, but it shows that my interpolation does something legit.

What you're seeing is acutance, as described in the wikipedia article on bicubic interpolation.

However, for downsampling an image, cubic interpolation is perhaps not the ideal choice, but makes more sense for upsampling.

Note that if you're downsampling by more than a facor of 4, there will be many pixels in the original image that are simply never used in the downsampled version.

The easiest way to get around this is to downsample by iteratively taking the mean value of 4 neighbouring pixels, until you're within a factor of 2 of the desired image size. Then you can do the last part by cubic interpolation.

If you're want high quality resampling, try taking a look at the ImageMagik documentation for inspiration. In particular they recommend using windowed filters for downsampling, such as Lanczos or Sinc.

I'm not sure what the problem is, but my first guess is that somehow
the 16 pixels that the code is pulling from the image is getting interpreted in the wrong order.

In other words, perhaps the block of code currently

val b_1 = q(x1, x2, x3)(extracted(0), extracted(1), extracted(2), extracted(3))
val b0  = q(x1, x2, x3)(extracted(4), extracted(5), extracted(6), extracted(7))
val b1  = q(x1, x2, x3)(extracted(8), extracted(9), extracted(10), extracted(11))
val b2  = q(x1, x2, x3)(extracted(12), extracted(13), extracted(14), extracted(15))

needs to be something else, perhaps something like

val b_1 = q(x1, x2, x3)(extracted(0), extracted(4), extracted(8), extracted(12))
val b0  = q(x1, x2, x3)(extracted(1), extracted(5), extracted(9), extracted(13))
val b1  = q(x1, x2, x3)(extracted(2), extracted(6), extracted(10), extracted(14))
val b2  = q(x1, x2, x3)(extracted(3), extracted(7), extracted(11), extracted(15))

?

That sort of shuffling pixels around would give results similar to what you are seeing -- the output looks fine in "smooth" regions, but the results are more obviously wrong in non-smooth regions.

This "nonscaling" test (run the "interpolation" to generate an image of exactly the same size) is a great first test -- it should give pixel-for-pixel exactly the same output image as the input image.
(Tests that expand or shrink the image, it's not as clear exactly what the output image ought to be).

I would also consider running a nonscaling test (or perhaps even a scaling test) on a very simple image, small enough that a human can single-step through every pixel in the image -- perhaps something like

0 0 0
0 5 6
0 8 9

Once you get this "nonscaling" test working, I agree with Robert that to really test the interpolation, you need to also run some tests that expand or shrink the image.


When to use the DTFT vs the DFT (and their inverses) in analysis? In many of my readings, whenever some author speaks about working in the frequency (transform) domain (of a digital signal), they often times take the DFT, or the DTFT, (and of course their corresponding inverses). Different authors will tend to work with one or the other.

I have not been able to really ascertain a particular pattern regarding this. In that, why would you pick the DTFT over the DFT or vice versa in explaining algorithms? Where does one help you over the other?

DFT and DTFT are obviously similar as they both generate the fourier spectrum of time-discrete signals. However, while the DTFT is defined to process an infinitely long signal (sum from -infinity to infinity), the DFT is defined to process a periodic signal (the periodic part being of finite length).

We know that the number of frequency bins in your spectrum is always equal to the number of samples processed, so this also gives a difference in the spectrums they produce: the DFT spectrum is discrete while the DTFT spectrum is continuous (but both are periodic with respect to the Nyquist frequency).

Since it is impossible to process an infinite number of samples the DTFT is of less importance for actual computational processing; it mainly exists for analytical purposes.

The DFT however, with its finite input vector length, is perfectly suitable for processing. The fact that the input signal is supposed to be an excerpt of a periodic signal however is disregarded most of the time: When you transform a DFT-spectrum back to the time-domain you will get the same signal of wich you calculated the spectrum in the first place.

So while it does not matter for the computations you should note that what you are seeing there is not the actual spectrum of your signal. It is the spectrum of a theoretical signal that you would get if you periodically repeated the input vector.

So I would assume in the literature you were mentioning, every time it is important that the spectrum you're working with is actually the spectrum and disregarding the computation side of things, the author would pick the DTFT.

The DFT is the transform of limited no of samples of a periodic signal.
the DTFT is a transform of the entire sample signal from $-\infty$ to $+\infty$ and input so not necessarily periodic.


How to determine the filter type of a discrete frequency spectrum I have the discrete impulse response of a filter and i want to determine which kind of filter it is by using the DFT.

The impulse response is:

h0 = [0,1,2,1,0]

which I zero-pad for a better result

N0 = 10;
h1 = [zeros(1,N0),h0, zeros(1, N0)];

Then i created the N-Point DFT of it using FFT

H1 = fft(h1);

The picture shows the plot of abs(H1)

Now i know that k1 = 0 represents the steady component and at first I thought that this is the spectrum of a band-stop filter but now that I read a little bit more about the DFT I'm not quite sure because the DFT is always mirrored for real signals.

So my question is: Do I have to use only one half of the DFT spectrum to determine the filter characteristics? (Which would mean that this is a low-pass filter?)

The $k_1$ indexes you get from the FFT correspond to frequencies 0 to $2\pi$ (actually, 1 sample short of $2\pi$). Note that frequencies in the range $\pi$ to $2\pi$ also correspond to the negative half of the specturm ($-\pi$ to 0).

You can get a nice visualisation with frequency 0 in the middle by applying fftshift to the result (Octave or Matlab assumed).

You can also check freqz, which plots a nice frequency and phase repsonse (try freqz([1 2 1], [1])). Or try plot(abs(freqz([1 2 1],[1]))) to get a linear plot rather than log.

What you are getting is a low-pass filter (gain 4 at low frequency, 0 at high frequency).

In VERY general terms: You can know beforehand it is a LPF because the impulse response is a single bump. A high-pass filter will have some sort of antisymetric shape (like [1 2 -2 -1]), and a band-pass or band-reject will have some oscillations in the impulse response.

You are getting a band pass filter. The coefficients of the filter are the ones that you are getting from the impulse response. The frequencies are folding around the f/2 or "the 0", where the imaginary values of the transform are anti symmetric (negative for half of the values).


Filtering on a time signal I am new to DSP, and am mostly self taught. I have a couple of questions, that may be too basic, but I couldn't really find an answer online.

I have a signal obtained by plotting the motion of a human in 2D from a fixed point of reference (origin).

The data is noisy, so I wanted to use a filter to get a smoother signal.
I want to filter out different motions. The intuition is that different kind of motions would occur at different frequencies.

From what I googled, and read online; it seems I need a bandpass filter like butterworth. However, after reading a few tutorials and book chapters, I have two questions:

My signal is a time-signal (displacement vs time), but the filter takes the input in terms of cut-off frequency. Does this mean I have to convert the signal into frequency domain, or does the filter take care of it? I am bit confused on how this actually works.
On a related note, I don't know the frequency of each motion that I am trying to filter out i.e. the intuition is that the frequency for different motions be different, but how do I determine by cut-off values? Is it usually done by trial and error, or is there a more scientific way to accomplish this?

My signal is a time-signal (displacement vs time), but the filter takes the input in terms of cut-off frequency. Does this mean I have to convert the signal into frequency domain, or does the filter take care of it? I am bit confused on how this actually works.

You design or specify the filter in the frequency domain, but the actual filtering can be directly accomplished in time.  There are FFT based filtering algorithms that are typically more efficient, but there is no requirement to transform your signal to the frequency domain.

On a related note, I don't know the frequency of each motion that I am trying to filter out i.e. the intuition is that the frequency for different motions be different, but how do I determine by cut-off values? Is it usually done by trial and error, or is there a more scientific way to accomplish this?

Anything that doesn't require trial and error is suspect, scientific methodology wise. I suggest that you try to be more specific about what filtering out different motions actually means.  Some kinds of phenomena are a superposition of features which are often separable, of which some can be "filtered out" and others don't lend themselves to that.

You might want to use a tool known as the Short Time Fourier Transform (STFT) and see if you can eye ball any patterns associated with desired/undesired motion.

Finally, try Google Scholar.

Few points to answer the questions

Editing the first answer from the comments of Marcus Muller

The transform to frequency domain and back always happens on parts of the signal, padded to specific lengths, and thus, the transform is part of the filtering algorithm itself, hence no need to explicitly transform time signal to frequency domain before filtering.
To find out the cutoff frequencies corresponding to each of the motions, you can see frequency spectrum of the input signal and try to figure out which frequency components are dominant during which type of motion. this will give you a fair idea of cutoffs for designing band pass filters.


What Is the Cause for Poor Results with Adaptive Thresholding? I try to run text cleaning process based on the camera captured image.

By using OpenCV adaptive threshold method with parameter (255.0D, AdaptiveThresholdType.MeanC, ThresholdType.Binary, 15, 10.0D), almost all gives satisfactory output but below example gives poor cleaning result (background is not removed well).

(image from http://www.fmwconcepts.com/imagemagick/textcleaner/index.php)

Question 1. I'd like to know which factor cause poor result in this example.

My Question 2 is, which approach can descriminate my attempt(OpenCV) is failed based on the above result? Or, which approach can help to decide whether AdaptiveThreshold will fail or not before running it?

IMO, the pass/fail factor can be a number of background particles but I don't know how to measure it.

I tried Erode method in OpenCV to remove particles, but it gave over whitened result and made too skinny text.

Some guidelines for using Thresholding:

Stretch the image to use the whole Dynamic Range (DR).
Apply some Denoising (Very very gentle).
Median with small radius would be a good idea.
Unless you hand tweak the Threshold, Otsu's Method generally yields good results for this kind of tasks (Text on background).
If one use Adaptive Local Methods (Mean / Gaussian based) the Radius is another hyper parameter to tweak.

Try using normal Thresholding function in OpenCV with type THRESH_BINARY.


Autocorrelation of discrete signal $\sin{2\pi fn}$ Determine the autocorrelation $r_{xx}[m]$ of the discrete signal

  $$x[n] = (\sin2\pi fn).$$

  where $n$ and $m$ are integers.

Using the definition I get

$$\begin{align}
r_{xx}[m] &= \sum_{n=-\infty}^{\infty}x[n]x[n-m] \\
&= \sum_{n=-\infty}^{\infty}\sin(2\pi fn) \sin\big(2\pi f(n-m)\big) \\
&= \sum_{n=-\infty}^{\infty}\sin(2\pi fn) \sin\big(2\pi fn - 2\pi fm\big) \\
\end{align}$$

but I can't seem to figure it out from here. I've tried using different trigonometric identities without result. I'm guessing it's something simple I'm missing.

well, there's an old trig identity that you learned pre-calculus:

$$ \sin(\alpha + \beta) = \sin(\alpha)\cos(\beta) + \cos(\alpha)\sin(\beta) $$

or, alternatively

$$ \sin(\alpha) \sin(\beta) = \tfrac12 \big(\cos(\alpha-\beta) - \cos(\alpha+\beta)\big) $$

and you will need to use one of those.

but you have a bigger problem because $x[n] = \sin(2\pi f n)$ is a finite power, infinite energy signal.  the definition for autocorrelation is different, i believe for power signals than for energy signals.

the equation for crosscorrelation for power signals is, i believe,

$$ r_{xy}[m] \ = \ \lim_{N \to \infty} \ \ \tfrac{1}{2N+1} \sum_{n=-N}^{N}x[n] \overline{y[n-m]} $$

for autocorrelation of a purely real signal

$$ r_{xx}[m] \ = \ \lim_{N \to \infty} \ \ \tfrac{1}{2N+1} \sum_{n=-N}^{N}x[n] x[n-m] $$

Correlation is the normalized form of covariance.  For zero meaned signals, a correlation of 1 means the two signals are proportional to each other, and -1 means they are negatively proportional.  A correlation of 0 means they are orthogonal, i.e. linearly independent.

From this, it is fairly easy to intuit what the answer of the OP's original continuous version of the question was.  (I don't know why it was changed it to the discrete case.)  Since $fl$, since converted to $fm$, is multiplied by $2\pi$ it is in units of whole cycles.  Therefore, when $fl$ is an integer, the autocorrelation is 1.  At half values, the autocorrelation is -1.  At the odd quarter values, the autocorrelation is 0.

So the answer ought to be:
$$ r_{xx}(l) = \cos(2\pi fl)  $$

If you use r b-j's second trig identity, and his second autocorrelation defintion, this result is fairly easily obtained.

I am sure the former is correct, but not sure about the latter.  According to Wikipedia's article the definition for autocorrelation is:

$$ \rho_{xx}(m) = \mathbb{E} \left[ ( X_n - \mu_x)( X_{n+m} - \mu_x) \right] / {\sigma_x}^2 $$

That's as far as I am taking this.

Hope it helps.

Ced


Implementing short-time Fourier transform I have a $512$ sample frame, lets call it X. The frame is further divided into four sub-frames of $128$ samples each. Each of these sub-frames have a different frequency.

For STFT, I use MATLAB's gausswin(N) function where N is equal to $128$ in my case to generate a Gaussian window, let's call it W. I zero-pad my windows W to make it equal to $512$ samples.

I now multiply my original frame X point by point with W to get y. I then take FFT of y to get Y. It works till here.

My question is, naturally the next logical step should be to shift my zero-padded window W to right by 128 samples and repeat the multiply and FFT procedure again.

But how do I accumulate these results?
Concretely, how do I convert it into a time-frequency plot?

admittedly, i have never used the buffer() function in matlab.  for doing STFT analysis and processing, i have just broken it into frames myself with code similar to this answer.  it's easier to see what's going on.

while you can use a gaussian window for analysis (i found it mathematically useful in this paper, email me if you want a pdf copy), you must use a complementary window, one that adds to 1 with the windows of its adjacent frames, for the reconstruction of the signal.  at least if you want this to pass the signal unchanged when the processing is set to "null".

To make things easier, I suggest you to use MATLAB function
buffer.
Regarding storage of your results I recommend to save it into matrix with dimensions  no_frames and $512$. So if your current frame index is i, then do:

S(:,i) = spec;

Where S is your output matrix and spec is the spectrum calculated for a given frame. Please mind that the S(:,i) will assign spec as i'th column in resulting matrix. If you need to transpose results, then just do S(i,:) or flip the spectrum by doing spec'.