Finding the frequencies given in a wav file I have a wav file containing a series of chords, each containing 3 or 4 notes, played concurrently.

I need to find the frequencies of each note in each chord, and I can use matlab to do any computation.

Currently I think I need to find the spectrogram of the wav file and convert the frequencies from normailsed frequencies to normal frequencies (I know the sampling rate). Here is the spectrogram:

Is this a good approach to solving this problem?

I've already tried it, but the spectrogram isn't particularly informative. Is there an alternative method that I'm ignoring?

Polyphonic pitch estimation is still a difficult problem.  Search on the MIREX forum for current research papers.  Simplistic spectrum based approaches usually require component pitches with extremely simple timbres (nearly pure unmodulated sine waves, etc.) rarely found in real-life polyphonic music.  More typical polyphonic music can contain a dense set of overlapping overtones plus other transient noises.  More complicated timbres creating polyphonic chords may require a statistical pattern matching based approach (against a database of potential instrument sounds), and/or machine learning algorithms, in addition to any signal processing for spectrum and envelope information.


How to Clean Black and White Images Which Are Results of Edge Detection I want to extract text from some images with messy backgrounds.

First step was switching to gray scale, applying an edge detection algorithm, and then cleaning up a little bit.

Here's an example of the result:

I want to use the tesseract OCR to extract the text from it, but it's not good enough to find it, as the background is quite messed up still. (I tried with a cleaner version and it's doing ok).

So, could anyone recommend an algorithm, or a type of algo, that would help me clean up this image and just get the text? (it will always have this font).

One idea I had was to try to remove all the clumps that are "larger" than the font, which is pretty thin.

ps: hopefully it's the right stackexchange site.

Turns out, the trick to solve this was just to realize that from frame to frame the background was changing, but the text was staying constant. So by doing an | operation on the pixel values of a few consecutive frames, it was crazy simple to get only the text.


Why does this video's audio change depending on the playing device? Yesterday I noticed that if I played this video from my smartphone on speaker, one could barely hear anything. On the other hand, if I used earphones I could hear it perfectly.

Then, I tried with several devices (other phones, my PC's speakers, a Bluetooth speaker) and on some of them one could listen to the audio with a really high quality, but on others it was very low and like 'echoish'.

I believe there should be a signal processing explanation to this phenomenon. I don't really understand what's going on, but I suspect it has something to do with stereo and mono audio. Sound quality was good on those devices that were stereo (speakers from PC, earphones, headsets) but on those that were mono (speaker of a cellphone, Bluetooth speaker) I could barely hear anything.

What is really going on? Is this phenomenon a mistake from the video uploader or could it be an intended effect?

The 2 stereo channel might be out of phase (near opposites, thus cancelling each other out when mixed to mono).

The human ear-brain doesn't just sum stereo left-right to the 2 ears, but instead uses any phase differences to help determine directionality.


Which antialiasing filter before equispaced sampling I am using MIT-BIH arrhythmical database where I have a digital signal of 1200 Hz so 1200 samples per second.
This means that analog filters have already applied to remove the frequencies over Nyquist frequency, so no aliasing.

However, I want to take equispaced sampling, every two sample, simply by the following in Matlab

data([0:2:1200]);

I am reading Andre Quinquils' book Digital Signal Processing Using Matlab 2008:

  It is always necessary to use an anti-aliasing filter before the
  sampling stage in order to avoid any spectral aliasing risk and to set
  an appropriate sampling frequency. In practice, a causal approximation
  of this ideal filter is used. Thus, depending on the chosen filter
  synthesis method, some imperfections are introduced, such as a
  passband amplitude ripple, a transition band and a stopband finite
  attenuation.

Does this mean that I need to apply a new anti-aliasing filter before the sampling stage in oder to avoid aliasing?
I think the sampling stage here is the equalspaced sampling.
I have not applied any new special anti-aliasing filter.

Which antialiasing filter can you use before equalspaced sampling stage?

Yes, you have to apply another filter before downsampling in order to avoid aliasing. Your original signal has been acquired with a sampling rate $f_\mathrm{s} = 1.2\,\text{kHz}$. Taking every second sample effectively means that you divide the sampling rate by two, so your new sampling rate $f_\mathrm{s}'= f_\mathrm{s}/2 = 600\, \text{Hz}$. To avoid aliasing you have to make sure that the signal contains no frequencies above $f_\mathrm{s}'/2 = 300\, \text{Hz}$ before you apply the downsampling. If it is certain that your signal does not contain any frequencies above $f_\mathrm{s}'/2$ then you don't need a filter.

The anti-aliasing filter can be implemented by an FIR filter. Its properties like stopband attenuation and width of transition band depend on your signal: if you have spectral components near to 300 Hz, the filter needs to be quite steep and it will consequently have a high order. If hardware implementation is an issue, have a look at Cascaded Integrator-comb filters.


getting different output when using FFTW3 Library than using Rosetta FFT I have the following array containing sine waveform that I want to compute FFT:
{ 0.0, 0.114124, 0.370557, 0.0, -0.576201, -0.370557, 0.0, 0.0 };

When using FFTW3 Library and running those 2 functions:
plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );
fftw_execute ( plan_forward );

I get the following output:
-0.462077
0.918921
-0.946758
0.233480
0.050790
0.000000
0.000000
0.000000

When using Rosetta implementation to FFT I get the following output:
-0.462077
0.918921
-0.946758
0.233480
0.050790
0.233480
-0.946758
0.918921

As you can see, some of the output is simillar between the 2s and some are different, I can't figure out why.
Here is my full code:
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
#include <complex>
#include <iostream>
#include <valarray>
using namespace std;

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
    const size_t N = x.size();
    if (N <= 1) return;

    // divide
    CArray even = x[std::slice(0, N/2, 2)];
    CArray  odd = x[std::slice(1, N/2, 2)];

    // conquer
    fft(even);
    fft(odd);

    // combine
    for (size_t k = 0; k < N/2; ++k)
    {
        Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
        x[k    ] = even[k] + t;
        x[k+N/2] = even[k] - t;
    }
}

int main(void)
{
    int i;
    double y;
    int N=8;
    double Fs=30000;//sampling frequency
    double  T=1/Fs;//sample time
    double f=5000;//frequency
    double *in;
    fftw_complex *out;
    double t[N-1];//time vector
    fftw_plan plan_forward;

    in = (double*) fftw_malloc(sizeof(double) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    for (int i=0; i< N;i++)
    {
        t[i]=i*T;
        in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
        double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window
        in[i] = multiplier * in[i];
    }

    Complex test[] = { 0.0, 0.114124, 0.370557, 0.0, -0.576201, -0.370557, 0.0, 0.0 };
    //const Complex test[] = (Complex*)malloc(sizeof(Complex) * N);
    //CArray data(test, 8);

    printf ( "\n" );
    printf ( "  Input Data:\n" );
    printf ( "\n" );

    for ( i = 0; i < N; i++ )
    {
        printf ( "  %4d  %12f\n", i, in[i] );
        test[i].real(in[i]);
    }
    CArray data(test, 8);

    fft(data);
    plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );
    fftw_execute ( plan_forward );

    printf ( "  output Data from FFTW3 Library:             output Data from my FFT implementation\n" );
    for ( i = 0; i < N; i++ )
    {
        printf ( "  %4d  %12f           %12f\n", i, *out[i], data[i].real());
    }
    return 0;
}

If the input is real valued than the output of the FFT is conjugate symmetric (Hermitian Symmetry). That means that the data in the second half of the output is redundant and the FFTW implementations takes advantage of that fact and only returns a complex vector of length $N/2+1$, i.e. 5 in your case. The Rosetta implementation just runs a standard complex FFT.
I suggest reading the documentation: http://www.fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html
Please note that your output is complex but you are only printing the real part.


Contour based blurring Suppose I have an image of a cat and a binary mask image corresponding to cat's location as shown

I would like to blur the image region surrounding the cat but as a function of distance from the cat's boundary pixels. Essentially, the further a pixel is away from the boundary, the greater the degree to which it gets blurred. The ideal distance measure would be iso-contours which extend radially from the cat's boundary pixels, but this is possibly complicated. I would settle for Euclidean distance for now.

Any idea how I can accomplish this ?

Compute the distance transform. You can do that very efficiently. This will give you the Euclidean distance to the closest boundary pixel, as well as the index of that pixel. Then, normalize all the distances in the mask so that the maximum becomes 1. This is required for the next stage: At each point initiate a Gaussian kernel. The kernel size $h$ (and similarly the standard deviation $\sigma$) should be a function of the weight, obtained from distance normalization. You can either use a coefficient - linear proportionality - or a complex kernel, depending on the blurring requirements.

Finally, blur the signal at each pixel using the generated mask. Each pixel gets a unique mask and thus a different blurring. Note that, some of the pixels could get the same blur.

You only need to adjust the kernel size $h$, as for typical Gaussian kernels, $\sigma$ is a function of $h$. The function is simply a lookup table:

   $h$  $\sigma$
    3    (0.65)
    5    (0.87)
    7    (1.43)
    9    (1.88)
    11   (2.31)


What is the optimal adaptive grid for calculating a DFT using a fixed number of sampling points? I'm currently facing the following problem: I want to approximate the Fourier transform $F(\omega)$ of a (let's say, $L^2(\mathbb R)$) function $f(x)$ by calculating the discrete Fourier transform, using only a fixed number of sampling points. These points are allowed to lie arbitrarily far away from the origin. In addition, we don't care about computational costs, so the grid can be chosen completely free.

What is the optimal (with respect to minimizing the error) choice for the grid?

I guess that this choice is dependent on the bandwidth $B$ of $f$. Does this mean that Nyquist-Shannon forces the grid to being broader than $2*B$?

Also, I found literature about sparse Fourier transform sampling in case of known sparsity in $F(\omega)$ (compressed sensing). But since my signal does not appear to have sparse support in the frequency domain, I don't see how these methods would be applicable in my case.

This can be split up in two questions for the sampling frequency and the window.

1) The required sampling frequency is given by the Nyquist theorem to be $f_s>2*f_{max}$ with $f_{max}$ being the largest frequency in your signal. When sampling actual signals with an ADC there is usually a low pass filter to ensure an upper frequency limit.

2) The function $f(x)$ is usually defined to $-\infty<x<\infty$. However, you can only pass a finite portion to the DFT. The DFT is not a time-discrete Fourier Transform but a Fourier Series. Hence, it is defined on periodic signals.

If your signal is actually periodic then it is favorable to perform the DFT over whole periods. From the DFT's point of view, the signal you hand to it is repeated infinitely and when you pass it an incomplete period, there will be border effects resulting in frequencies showing up that are not in the original signal.

If your signal is time limited, all of it should be contained in the vector you pass to the DFT.

If the signal is neither time limited nor periodic you would theoretically need a DFT with infinite length, simular to the actual Fourier Transform which is defined for $-\infty<x<\infty$. Since this is not possible, Spectral Density Estimation techniques can be employed. The data can then be visualized with Spectogram plots showing the frequency content over time. It depends on what you are trying to achieve, really.


What is the estimator used in this paper for system identification? The paper System Identification using Symbolic Chaotic
Sequence proposes EM-UKS estimator for system identification of a linear FIR channel when excited by non linear input.

In Fig 3 of the paper there is a comparative plot of system identification of FIR model

$$ y[n] = H^T x[n] + w[n] $$

Where $x[n]$ is a chaotic input and $w[n]$ is AWGN with unknown variance.

Performance is compared using CRLB vs MSE when the FIR system is driven by (a) chaotic deterministic numerical input, (b) chaotic quantized input (called as symbolic input), (c) filtered Gaussian input, and (d) PRBS input.

Q1: What is the estimation technique when the input is white Gaussian and PRBS ? Is it Ordinary Least Square as they mention that they are comparing with the CRLB of non-blind. So, OLS is a non-blind technique.

This is how I implemented when the input is white noise:

 h = [1    0.65   -0.2];  % MA channel coefficients
    x = randn(1000,2);  % white noise input
    u = filter(h,1,x);
    y = awgn(u,10,'measured'); % adding noise w of SNR =10dB
    X = [y,ones(n,1)];
    [b,C]=lscov(X,y)

Q2: I am not sure how to implement using PRBS. Can somebody please provide the code?

Q3: Is the system identification using EM-UKS when the input is chaotic and symbolic proposed by the Authors blind or semi-blind and why?

Interesting paper!

Q1: yes, seems like OLS, non-blind. Usually OLS is non-blind. Some sort of regulated least squares can be used in a bayesian ML-setting, blind. This paper looks interesting: http://www.cs.berkeley.edu/~jordan/papers/lindsten-etal-sysid12.pdf

Q2: PBRS: https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence. There is a lot of literature on this, e.g. Coleman Brosilow, Babu Joseph, "Techniques of Model-Based Control"

Q3:they say it is semi-blind since the dynamics of the sequence generator is known.


Discussion - Integrating Accelerometer data to Position data from Frequency domain FFT I came across one blog who suggested to convert from accelerometer data to position data from frequency domain. Here are the steps that the original individual suggested:

  Remove the mean from the accel. data
  Take the Fourier transform (FFT) of the accel. data.
  Convert the transformed accel. data to displacement data by dividing each element by -omega^2, where omega is the frequency band.
  Now take the inverse FFT to get back to the time-domain and scale your result.

I am wondering if anyone have ever used this method before, and may I know your algorithm? I am not particularly sure on what does it means by Step 3 - frequency band?

p.s. I hope this question is relevant here.

The point is that acceleration is the second derivative of displacement. So to get the displacement, one should integrate acceleration twice. To do this, one way is to numerically integrate. A more interesting approach is to use the time-integration property of the Fourier transform:
$$\mathcal{F}\left(\int_{-\infty}^t x(\tau)d\tau\right)=\frac{X(j\omega)}{j\omega}$$
Note that we require $X(0)=0$. So division by $-\omega^2$ means twice integration in time.
Assuming the sampling frequency is fs and a is accel. data, you can use something similar to the following code:

 l = length(a); % assume l even
 w = 2*pi*fs/l*(-l/2:l/2-1);
 A = fftshift(fft(a));
 D = A./-w.^2;
 d = ifft(ifftshift(D),'symmetric');


Algorithm/Code for finding surface shape of an image? Is it possible to find the surface shape of an image as mentioned in the reference question ?
If I find the unknowns of the equation , then what should be my next step in determining the shape of the surface/image.
Kindly mention the algorithm/Code or the procedure required for this purpose.

Thankyou.

The 3x3 linear system appears to be that for photometric stereo, where the $i$-th row is $$l_i\cdot n = b_i,$$
with

$l_i$ is the 3D light source vector, e.g., $l_1$ the first row of the matrix, (0.2425, 0, -0.9701)
$n = (x_1, x_2, x_3)$ is the unknown surface normal vector, or "shape"
$b_i$ is the brightness observed when the surface was illuminated with $l_i$

Solving this 3x3 linear system gives the normal $n$ for a single surface point ("middle" according to "Example 1").  (You may wish to divide $n$ by its length $\sqrt{x_1^2+x_2^2+x_3^2}$ so that it is unit length.)

Usually, by surface shape in this context of photometric stereo one means the height function $z(x,y)$, where we associate a height to each pixel $(x,y)$.  To obtain this, we first need to solve a related 3x3 system at each pixel location $(x,y)$ (not just the "middle"), so that we have a surface normal $n(x,y)$ everywhere.  Then we solve for that function $z$ whose normal function is that $n$ we obtained via many 3x3 solutions above.  Finding $z$ requires we impose integrability on vector field $n$ (because modeling error and noise may mean there exists no $z$ having normals $n$).  However, I am digressing from your original "Example 1", which seeks a single surface normal.

See also: http://pages.cs.wisc.edu/~lizhang/courses/cs766-2008f/syllabus/10-09-shading/shading.pdf and https://en.wikipedia.org/wiki/Photometric_stereo

Note that the above assumes diffuse reflection (no surface glossiness), surface smoothness, and no shadows nor occlusion.  Light sources must be distant if we keep them constant over image.


GCC-PHAT always peaks at zero when estimating delay offset My question is that when I used gcc-phat to estimate delay between two audio signals, it always returns zero offset. However normal cross correlation Xcorr gives a good result. I don't know the reason.
It is the same case as below two:
GCC-PHAT (Generalized cross correlation MATLAB)
GCC-PHAT (Generalized cross correlation) always peak at delay=0 on real audio signal

I have tried to remove the DC part of signals but it doesn't work.
Any suggestion is appreciated!

This code seems to work (included also below).
Pythonizing the first example link you include seems to do the right thing with that code.
Fs = 8000
dt = 1/Fs #0.125e-3
f1 = 100
tdelay = 0.625e-03 # try different values
t3 = np.linspace(0.0,1.0,Fs)
x3 = np.cos(2*np.pi*f1*t3);
x4 = np.cos(2*np.pi*f1*(t3-tdelay));
gcc_phat(x4,x3)

yields:

(0.375,  array([-0.02514735, -0.02196879, -0.01847165, ..., -0.0303065
,
-0.02794452, -0.02514735]))

(which gives one minus the actual delay).
After playing with your wav file (thanks!), I think the problem is that your recording equipment has too much correlated noise across all four channels.  That's giving a GCC-PHAT peak at zero.
If I try to drown out the noise by adding even more (uncorrelated) noise, then I get a better result: -5 samples vs -4 samples from XCORR.

The four plots are:

The section of one channel that I'm looking at.

The GCC-PHAT central part when uncorrelated Gaussian noise is added to the data.

The GCC-PHAT central part when no noise is added to the data.

The XCORR result of the noiseless data.

Code to do this is below and on GitHub.  I originally though it might be possible to just filter the data, but that doesn't help at all.
import numpy as np
import statistics
from scipy.io import wavfile
import matplotlib.pyplot as plt
from IPython.display import Audio
from scipy.signal import kaiserord, lfilter, firwin, freqz

samplerate, data = wavfile.read("Q69905.wav",'rb')
samples = np.arange(20000,30000)

fir_filter = firwin(123,0.75)
channel_1 = lfilter(fir_filter,1, data[samples,0])
channel_2 = lfilter(fir_filter,1, data[samples,1])
channel_3 = lfilter(fir_filter,1, data[samples,2])
channel_4 = lfilter(fir_filter,1, data[samples,3])

noise_1 = np.random.normal(0,1000,len(channel_1))
noise_2 = np.random.normal(0,1000,len(channel_2))
noise_3 = np.random.normal(0,1000,len(channel_3))
noise_4 = np.random.normal(0,1000,len(channel_4))

print([statistics.mean(data[:,0]),  statistics.mean(data[:,1]), statistics.mean(data[:,2]), statistics.mean(data[:,3])])

# delay, gcc = gcc_phat(data[samples,0].astype(float)+10, data[samples,2].astype(float)+12, interp=1)
delay, gcc = gcc_phat(channel_1 + noise_1, channel_3 + noise_3, interp=1)
delay_no_noise, gcc_no_noise = gcc_phat(channel_1 , channel_3 , interp=1)

plt.figure(figsize=(20,30))
plt.subplot(4, 1, 1)
plt.plot(data[samples,0])
plt.subplot(4, 1, 2)
plt.plot(np.arange(-10,10),gcc[9990:10010],'.') # [9950:10050]
plt.subplot(4, 1, 3)
plt.plot(np.arange(-10,10),gcc_no_noise[9990:10010],'.') # [9950:10050]
plt.subplot(4, 1, 4)
lags,c, line, b = plt.xcorr(channel_1,channel_3)
plt.plot(lags,c,color='r')

print('GCC-PHAT: ' + str(delay))
print('XCORR: ' +  str(lags[np.argmax(c)]))
Audio(channel_1 + noise_1, rate=44100)

Code From GitHub
"""
 Estimate time delay using GCC-PHAT
 Copyright (c) 2017 Yihui Xiong
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
     http://www.apache.org/licenses/LICENSE-2.0
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
"""

import numpy as np

def gcc_phat(sig, refsig, fs=1, max_tau=None, interp=16):
    '''
    This function computes the offset between the signal sig and the reference signal refsig
    using the Generalized Cross Correlation - Phase Transform (GCC-PHAT)method.
    '''

    # make sure the length for the FFT is larger or equal than len(sig) + len(refsig)
    n = sig.shape[0] + refsig.shape[0]

    # Generalized Cross Correlation Phase Transform
    SIG = np.fft.rfft(sig, n=n)
    REFSIG = np.fft.rfft(refsig, n=n)
    R = SIG * np.conj(REFSIG)

    cc = np.fft.irfft(R / np.abs(R), n=(interp * n))

    max_shift = int(interp * n / 2)
    if max_tau:
        max_shift = np.minimum(int(interp * fs * max_tau), max_shift)

    cc = np.concatenate((cc[-max_shift:], cc[:max_shift+1]))

    # find max cross correlation index
    shift = np.argmax(np.abs(cc)) - max_shift

    tau = shift / float(interp * fs)

    return tau, cc

def main():

    refsig = np.linspace(1, 10, 10)

    for i in range(0, 10):
        sig = np.concatenate((np.linspace(0, 0, i), refsig, np.linspace(0, 0, 10 - i)))
        offset, _ = gcc_phat(sig, refsig)
        print(offset)

if __name__ == "__main__":
    main()


Spikes in time series I would like a simplified measure to quantify spikes in my time series, I have one series that has many and I would like to compare it to another that has very few. I have thought through several options for example fitting a sine wave and looking at the R values as a measure of smoothness etc but I thought that a measure of the Coefficient of variation of the differences in the series may be suitable. I was wondering if this was a good approach and also as the mean values of the differences are very small I was wondering if theres a common alternative for the CV which of course would become meaningless.

Also just to mention that there is no proportionality between the means and standard deviations of the data.

Another way to ask this question would be if anyone had suggestions for the inverse, a successful measurement of smoothness? or a way to summarise jerking in the data

I would try a median filter.

Let your original signal be $f[n]$.
Median filter $f[n]$ using $N$ pixels, where $N > 2 \times S + 1$, where $S$ is the maximum number of samples in the spike. The resulting signal, lets call it $g[n]$ should have all the spikes removed.
Find the absolute of the difference between the two signals, $h[n] = |f[n] - g[n]|$. This signal represents the spikes.
Count the number of positive transitions in $h[n]$ that are above a threshold. This is the number of spikes.

Example:

Original signal $f[n]$

Median filtered signal $g[n]$

Absolute difference $h[n]$


References and resources on hyperspectral (multispectral) image processing Due to my work, I just started studying hyperspectral or multispectral image processing.

I think I have enough programming background and studied mathematics for my degree - it was long time ago though. Since early this year, I have been exposed to image processing in my development work for manufacturing industry. And I got involved in a research project dealing with the subject above.

I wonder if anyone can recommend good references or resources for the subject. Especially, I would like to be familiar with that kind of image data. Any advice would be appreciated.

What kind of information are you after.  Are you looking for information on how to process the  images from a data processing point of view?

Or from a theory point of view?

Are you looking at satellite images?  if so land or water?

Proprietary software look at ENVI

http://www.exelisvis.com/language/en-us/productsservices/envi.aspx

for free software look at SeaDAS

http://oceancolor.gsfc.nasa.gov/seadas/


How to recognize speech bubbles in comic strips (Ideally OpenCV) I need to recognize speech bubbles to send to OCR and also to know where they are in the strip.
How could I do this with OpenCV?

I'm using OpenCV 2 and Python 3.6

As said in the comments an efficient way is to first detect letters, words and text with OCR. Then try to expand each text zone to its corresponding text bubble.

Depending on the text bubble design there are different approaches. However, a solution that could work well and be robust would be to perform edge detection on the near surrounding of the previously detected text zones and then keeping only the (radially) closest edges from the text zones(that should be the bubble borders).

This is just a basic idea, you may have to build a more sophisticated algorithm obviously!


Converting from RGB to LAB Color Space to Reduce Colors In many papers I see the following:

  First we have to convert image from RGB to LAB color space to reduce colors then algorithm can run faster.

I can't understand this. How can this happen? What I read is that Lab color space has support for more colors than the other. But I see people says converting to Lab reduces the colors to better calculation. Why and how does this happen?

You should read about LAB Color Space.

Its main advantage on RGB is by having 1 channel dedicated to the Luminosity of the Image and 2 others dedicated to the color information.

This means if you need to process only Luminosity / Color you can do that using LAB without having any effect on Color / Luminosity.


Framing an audio signal The first step  in MFCC is to split the signal into frames and  referring to this discussion on MFCC at this link (Help calculating / understanding the MFCCs: Mel-Frequency Cepstrum Coefficients) @pichenettles said that  the frames are usually  overlapping.
my question is why they are overlapping ?. And if they don't , for example
"frame 0 is from 0->440 and frame1 is from 440 ->880 " does this affect the process of feature extraction using MFCC method?

And also why a widowing function must be applied to the frame ?

Thanks.

The MFCC summary you link seems to leave out the typical windowing function applied before each FFT.  Segmenting longer data into shorter finite length FFT inputs does an implicit rectangular windowing, which causes the energy of the frequency of any signal not exactly integer periodic in the FFT length to be "spattered" into other distant frequencies of the FFT result.  A non-rectangular window function is applied to help reduce these rectangular window artifacts (sometimes called spectral "leakage").

A typical window function, such as Von Hann or Hamming, can be quite lossy at both edges of each window.  Overlapping windows adds a window (or more) centered closer to where other windows are lossy, thus helping preserve information about the original signal that would otherwise be lost or degraded.

Another reason is for time resolution.  The length of each FFT window might be chosen to give sufficient frequency resolution, however, due to the time-space resolution trade-off of FFT analysis, the resulting time resolution may not be adequate, and thus miss fast audio transients not being centered or isolated in one window.  Shorter FFT windows may not provide sufficient frequency resolution.  A greater density of FFT windows (using overlap instead of shorter windows) improves the capture and isolation of these shorter time domain events.  This greater density (more windows per unit time) can be accomplished by overlapping the windows.

If you don't overlap, then events near or at the boundaries will be severely degraded, and the probability of transients (such a consonants) being between windows, or not being isolated in a single window, is increased, thus reducing how well the MFCC can categorize different inputs.


Making audio clips comparable in terms of RMS I have a challenging question on how to normalize several audio clips so that they are comparable (similar) in RMS but without the clips clipping (having values over -1 and 1). Is there a systematic way to do this? One scaling factor must be found for all clips so that their L/R channels keep their ratio while at the same time making the different clips' RMS equal. Thanks in advance!

I wrote this simple Matlab function a while back, it operates on dBFS:

function y = change_rms(x, newRMS)

    currentRMS = mean(20*log10(rms(x)));

    d = newRMS - currentRMS;
    y = x * 10^(d/20);
end

You would call a function like this with the signal whose RMS you want to change as the first argument (consecutive samples along the first dimension, channels along the second) and the new RMS value in dBFS as the second argument. After processing all signals so that they have an RMS equal to your reference RMS, you check which signal contains the highest absolute value. Then devide all signals by that value. They will then all have the same RMS (i.e. mean RMS) and their internal level ratio between channels will still be what it was before changing their RMS.

Note however that the final RMS that all signals will have after this procedure will not be equal to the reference RMS you used when calling the function. This is of course due to dividing all signals by the maximum absolute value of the entire data set after using the function. Consequently, the final RMS will be determined by the signal that contained the maximum absolute value across the entire data set.

Here is the step-by-step guide:

Process all your audio files with the script provided above. Use the
same newRMS for all the audio files. newRMS can be an arbitrary
value.
Find the highest absolute value in the entire set of audio files.
Divide all audio files by that value.

After that, all audio files will have the same RMS. All audio files will have values in the [-1 1] range.


Pose estimation of colored checkerboard cubes I want to find the orientation of 2x2 checkerboard cubes, like these:

That is to say, I want to recover their positions and rotations within the image relative to the camera.

In my particular case, I care more about the rotation than the position (and the depth barely matters at all). I know the colors ahead of time, and can change them to be more convenient. I want to pick up the cubes, so a hand may be in the frame. The lighting may not be uniform, but the background should be simple (e.g. carpet, table, pavement).

I do have a partially working solution, which I've included below as an answer (it makes the question too specific to include it here). But I'm looking for suggestions on how to approach the problem and what techniques to apply.

Other things I've been considering:

Use a line detection algorithm to detect edges, and use those to seed hypotheses (instead of the center saddle point).
Do a cartesian-to-polar transform at the detected saddle points to turn the angled lines into straight lines (expensive?).
Score hypotheses by using them to unskew a face then dot product-ing that against a reference image of the face (but what about lightning changes?).
Random combinations of dilation and erosion throughout (invariably doesn't seem to help).

This is the partial solution alluded to in the question. Should I be using a totally different approach (probably)? Are there standard operations that do what I've cobbled together (definitely)?

Sample circularly at a radius of ~6 around each pixel and compute the sum of the absolute differences to the opposite side. This seems to do a pretty good job of picking up edges, and smoothing noise:

Do the same thing again but, instead of comparing across to the 180 degree counterpart, compute differences to the 90 degree counterparts. This still picks up edges, but especially picks up the saddle points at the centers of each face (less so when tilted, though):

Subtract the previous two transformations, cancelling the edges and keeping the saddle points. This seems to work quite well:

Find isolated local maxima in the saddle point image. At each one over some threshold, sample every pixel a distance of ~10 away. Find the alternating colors via a kmeans with k=2 on those pixels' colors. Use the resulting labels and a simple maximization algorithm to estimate the angles where the colors switch. If the kmeans match is terrible or there's lots of mismatches between the angles, discard the point as a false-positive::

For each remaining saddle point, offset by ~15 along each of its diagonal directions and do a flood fill against a finer version of the half-turn-edges image, made by sampling at a radius of ~2 instead of ~6:

For each flooded area, compute its contour and fit a rectangle. If the rectangle covers significantly more area than the contour, the flood fill probably escaped the face. For the flood fills that worked, use the furthest points along the diagonals (the support points) to estimate the position of the corners and the side-centers. If the two flood fills adjacent to a side succeed, average the points. When 2/4 or 1/4 of the flood fills for a saddle point failed, use the existing points to infer the missing ones:

Use the four estimated corners to compute a perspective warp. Try to refine the warp by tweaking the points and preferring lots of brightness near the sides and equators of the half-turn-difference transformed image... and otherwise lots of darkness in the interior.

(todo) Fuse the found faces into cubes. Recover the entire 3d pose.

From my experimenting, the weaknesses of this method are mainly:

Bad at finding saddle points on faces titled more than ~45 degrees.
Bad at refining the faces, unless the flood fill gets it extremely close to correct.
Bad at excluding close false positives, like if the saddle point on an edge is assumed to be a center.
Any blurring from motion ruins the saddle point discovery.
Doesn't infer poses from the faces, yet.

Also its kind of slow, running at about 1fps, but I haven't put much effort into optimizing it.

Any suggestions?


Factor $|a|^{-1/2}$ in definition of mother wavelets The mother wavelet is defined as

$$\Psi_{a,\tau}(t) = \lvert a\rvert^{-1/2}\Psi\left(\frac{t-\tau}{a}\right)\tag{1}$$

in which the factor

$$ \lvert a\rvert^{-1/2}\tag{2}$$

is supposed to serve that

$$\lVert\Psi_{a,\tau}\rVert = \lVert\Psi\rVert\tag{3}$$

is valid.

Does someone have an idea how $(3)$ is to be understood?

My answer is for real scale $a$ and the fact that wavelet transform is usually defined in $L_2$ with norm

$$||\Psi(\tau)|| = \int_\mathbb{R} \Psi(\tau)\Psi^*(\tau)\mathrm{d}\tau $$

So

$$||\Psi_{a,t}(\tau)|| = \int_\mathbb{R} \frac{1}{|a|}\Psi(\frac{\tau-t}{a})\Psi^*(\frac{\tau-t}{a})\mathrm{d}\tau$$

Set $\tau' = \frac{\tau-t}{a} \implies d\tau' = d\tau / a \implies d\tau = a \times d\tau'$

$$||\Psi_{a,t}(\tau)|| = \int_\mathbb{R} \frac{1}{|a|}\Psi(\tau')\Psi^*(\tau') \times a \times \mathrm{d}\tau' = \int_\mathbb{R} \Psi(\tau')\Psi^*(\tau') \mathrm{d}\tau' = ||\Psi(\tau)||$$

Update for the question "How does it not change the norm considering that it could introduce a negative sign with the norm being sign-less?"

We develop stuff for the case that scaling factor $a < 0$.

with $A > 0$
$$||\Psi_{a,t}(\tau)|| = \lim_{A \to \infty}\int_{(-A-t)/a}^{(A-t)/a} \frac{1}{|a|}\Psi(\tau')\Psi^*(\tau') \times (-|a|) \times \mathrm{d}\tau' \\
= \lim_{A \to \infty}\int_{(A-t)/a}^{(-A-t)/a} \Psi(\tau')\Psi^*(\tau') \times \mathrm{d}\tau' \\
= \lim_{A \to \infty}\int_{(-A+t)/|a|}^{(A+t)/|a|} \Psi(\tau')\Psi^*(\tau') \times \mathrm{d}\tau' = \int_{-\infty}^{+\infty} \Psi(\tau')\Psi^*(\tau') \times \mathrm{d}\tau'$$


Converting blobs to circles I have this binary image with some blobs of random sizes.

I would like to convert the blobs into circles of uniform diameter. it doesn't really matter what diameter but they need to be circles and uniform.
Here is what I am looking for.

I am working with tiff sequences in ImageJ. I tried erosion and dilation. I tried Ultimate point to reduce to a point and then dilate up but that creates a square. I tried distance map but didn't get very far. Although i am pretty sure distance map is the way to go. None of my current approaches are satisfactory.

If I could do

   #pseudo-code
   if area(blob) > x
   then make circle(diameter=y)
   else remove blob

would be even better.

If anyone has a solution using ImageJ or Python SciPy or even R, please let me know. I would very much appreciate that. Thanks.

Here is the procedure to follow:

Threshold your image
Obtain connected components - each distinct blob in the image gives you a component.
For each component, find the boundary pixels. Each boundary is represented as a set of points.
Fit one circle per each boundary point set. You could use a simple least squares fitting.

This should give you pretty good circles. You could do this in many languages, using the standard image processing libraries.


What does "kernel based" mean? In my thesis I try to explain what kernel based methods are, especially the meaning for object detection.

I know kernel based methods like Mean- and CamShift and I know how to use them. I understand how the shift work. But: What does the kernel do, what does he describe?

I know, wikipedia have articles about kernels but I still don´t get it. :(

Q1: What could be an subset for an image?

Q2: How does the kernel project the points in an image?

Q3: Could you give me an simple example to understand kernels?

Thank you in advance!

In general, a kernel is a function that acts as a parameter to some algorithm.

Filtering: For example, it's possible to call the impulse response of a filter $h[n]$ a kernel, so that it is the parameter that defines the filter operation:
$$
y[n] = h[n] * x[n].
$$

The use of the term kernel in the filtering context is much more common in 2D filtering or image processing.  The link talks about the kernel being a matrix, but really it's just a sampling of the function that is the "true" kernel.

PDF Estimation: Kernel-based methods are often used in other contexts, too.  For example, when estimating the probability density function of a random variable, kernel-based estimators are often preferable to simple histogramming.  In that context, there are many different possible kernels.

Machine Learning: Finally, another context for kernel based algorithms is in machine learning.  Here, we are interested in classification of an input into one of (possibly) many classes.  Again, the kernel is a function $k(\mathbf{x}_i,\mathbf{x}')$ that parametrizes the algorithm and there are many possible selections.


Distribution of the output of AWGN channel $$\mathbf{Y}=\mathbf{X}+\mathbf{W}$$

$\mathbf{X}$ is input, an $n$-dimension equiprobable vector in a finite subset of complex space $\mathcal{A}^n \subset \mathbb{C}^n$ with norm $\sqrt{\rho n}$ (constant power), e.g. BPSK, QPSK, 8-PSK, etc.

$\mathbf{W}$ is circularly complex Gaussian vector $\mathcal{CN}(0, \mathbf{I}_n)$.

I want to know the closed form of $P_{\mathbf{Y}}(\mathbf{y})$. And if possible, the $P_{\mathbf{Y}}(\mathbf{y})$ induced by the QAM modulation of $\mathbf{X}$ (identical average power).

I can numerically calculate
$$P_{\mathbf{Y}}(\mathbf{y}) = \sum_\mathbf{x}P_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})P_{\mathbf{X}}(\mathbf{x}) = \sum_\mathbf{x}P_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})\frac{1}{|\mathcal{A}|^n} = \frac{1}{|\mathcal{A}|^n} \sum_{\mathbf{x} \in \mathcal{A}^n} \frac{1}{\pi^n}e^{-\left\Vert\mathbf{y-x}\right\Vert ^2} \tag{1}$$

where $|\mathcal{A}|$ is the cardinality of set $\mathcal{A}$. But I want to avoid $\sum_{\mathbf{x} \in \mathcal{A}^n}$ because it is cumbersome.

If $\mathcal{A} = \mathbb{C}$ I can use the spherical symmetry trick that
$$P_{\mathbf{Y}}(\mathbf{y}) = \frac{1}{S} \times P_{\left\Vert\mathbf{Y}\right\Vert^2}(\left\Vert\mathbf{y}\right\Vert^2) \tag{2}$$

where $S$ is the surface of n-dimension complex sphere with radius $\sqrt{\rho n}$ and $P_{\left\Vert\mathbf{Y}\right\Vert^2}(\left\Vert\mathbf{y}\right\Vert^2)$ is noncentral $\chi$-squared distribution with $n$ degrees of freedom and noncentrality parameter $n \rho$. Thus I have a closed form of the distribution of $\mathbf{Y}$.

I think this must be a classical problem, it would be great if someone gives me a hint or points out a paper.

I see this question is open for a long time, I am not sure if I understand it correctly.
PDF of the symbols
The first probability is a $P_{\mathbf{Y}}: \mathbb{C}^n \to \mathbb{R}$, you are not considering phase noise, the probability  of the $n$ symbols are identically distributed, let's say with $P_{\mathbf{Y}_i}: \mathbb{C}^n \to \mathbb{R}$, since the joint probability of independent variables are the product of the individual variables this distribution can be written as
$$P_{\mathbf{Y}}(y) = \prod_{i=1}^{N}P_{\mathbf{Y}_i}(y_i)$$
Where $P_{\mathbf{Y}_i}(y_i) = \sum_{x\in \mathbf{X}_i}p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}/\pi$ which leads to
$$P_{\mathbf{Y}}(y) = \frac{1}{\pi^N}\prod_{i=1}^{N} \sum_{x\in \mathbf{X}_i} p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}$$
The term $p_{\mathbb{X}}(x_i)$ can be used to model probabilistic shaping, but for the original quesiton it will probably be a the inverse of the number of symbols in the constellation $1/c = 1/|\mathbf{X}_i|$, that by it's turn is the same for all symbols, thus
$$P_{\mathbf{Y}}(y) = \frac{1}{c^N\pi^N}\prod_{i=1}^{N} \sum_{x\in \mathbf{X}_i} p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}$$
The complex of evaluating this will be $c N$ exponentials. The number of additions and multiplications will be linear in $cN$ as well, differently from the initial formula that would compute $c^N$ terms.
For some applications you may want the logarithm of that expression, in such situations it is preferable to use the sum of logarithms
$$log\left(P_{\mathbf{Y}}(y)\right) = \sum_{i=1}^{N}\log\left(\frac{1}{c \pi} \sum_{x\in \mathbf{X}_i} p_{\mathbf{X}}(x_i) e^{-||y_i - x_i||^2}\right)$$
All of this is valid for QAM modulation as well, and don't have the exponential number of terms in the initial formula.
PDF of the energy
In the second part you mention the PDF of the energy of the symbols
That is the probability distribution of
$$ \sum_{i=0}^{N} \sum_{x_i \in \mathbf{X}_i} p_{\mathbf{X}_i}(x_i) ||x_i + w_i||^2$$
Assuming $w$ having a variance $\sigma$ for both real and imaginary parts, since the norm of a complex number $||a + 1j b||^2 = a^2 + b^2$, assuming uniform probability for the symbols in the constelation, we can divide the whole expression the symbols $\sigma$ (since the the formula is given for a covariance matrix being the identity). And the distribution of the energy will be of the form $\sum_{k=1}^{2N} X_{i}^2$ as assumed for the noncentral $\chi^2$ distribution with $2N$ degrees of freedom.
The noncentrality parameter given by
$$\lambda =\frac{N}{c \sigma^2} \sum_{x \in \mathbf{X}} ||x||^2$$
The number of degrees of freedom $k=2N$, because we separated real and imaginary parts.
Let $z$ be the energy of one sample, taking in to account that the $\chi^2$ distribution was normalized by dividing by the variance, the PDF of $z$ will be
$$p_{||\mathbf{Y}||^2}(z) =
\frac{1}{2\sigma} e^{\left(\frac{z}{2\sigma} + \frac{\lambda}{2}\right)}
\left(\frac{z}{\lambda \sigma}\right)^{(N-1)/2}
I_{N-1}\left( \sqrt{\frac{\lambda z}{\sigma}}\right)
$$
Note:
In spite of being a closed formula it depends on the computation of $I_{\alpha}(x)$, a modified Bessel function, that for positive integer $\alpha$ and positive $x$ reduces to a series with rational coefficients in $x$, thus
$$ I_{N-1}\left( \sqrt{\frac{\lambda z}{\sigma}}\right) = \sum_{m=0}^{\infty} \frac{\left(\sqrt{\frac{\lambda z}{\sigma}}\right)^{2m+N-1}}{m!(m+N)!} $$


How DFT is applied to Signals? I know the size of DFT is the number of samples used to do the calculation.
But the audio signal (and other signals) normally have many samples.
So, what does really happen, if I choose $2$ seconds of audio, which will have $2 \times 44100$ samples, and $1024$ as the size of DFT ?

First, the result of the DFT is the amplitude and phase of the frequencies which are stepped by $44110 / 2 / (1024 / 2) \approx 43\textrm{Hz}$.
So we get $0\textrm{Hz, } 43\textrm{Hz, } 86\textrm{Hz, } .... \textrm{to } 22050\textrm{Hz}$.
Is this correct?

Second, as I analysis the $2$-second signal, the calculation is done by using every $1024$ samples to calculate $512$ DFT bins, then what ? We can get $44100 * 2 / 1024 = 86$ results of DFT, we calculate the mean of every result to build the final result ?

In order to understand the effect of DFT length $N$ on the computed spectrum samples, you shall be aware of two things:

1- The interrelation between theoretical DTFT and practical DFT.

2- The windowing effect on the analysed samples

For the first relation, consider a causal, discrete-time signal $x[n]$ of finite length $L$, whose DTFT, $X(e^{j\omega})$ is:

$$  X(e^{j\omega}) = \sum_{n=0}^{L-1} x[n] e^{-j\omega n}$$

This is mainly a theoretical tool, it's a complex valued continuous-frequency function which can not be represented in a digital computer system which can only work with sampled data of finite bit length.

In order to solve this issue, an alternative definition of Fourier transform for discrete time signals have been proposed and called as the DFT (discrete Fourier transform) $X[k]$ given as:

$$  X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N} n k} ~~~,~~~ \text{for} ~~~ k=0,1,...,N-1$$

Where $N$ is the DFT-size and not necessaily the signal length (usually $N \geq L$ is taken). When $N=L$ it can be seen from visual inspection that:

$$ X[k] = X(e^{j \omega})|_{w = \frac{2\pi}{N} k}  =  X(e^{-j\frac{2\pi}{N} k})$$

That's to say, N-point DFT values $X[k]$ of finite length signal $x[n]$ are the uniform samples of the DTFT $X(e^{j\omega})$ of the signal $x[n]$ itself. Hence we maintain the view that DFT values are the samples of the DTFT evaluated at $w_k = \frac{2\pi}{N} k$ frequencies.

Now it's obvious that the more samples you evaluate (larger the $N$) on DTFT $X(e^{j\omega})$ the smoother representation $X[k]$ will have. In addition, the minimum number of DFT samples that could exactly represent DTFT is given by $N_{min} = L$. Furthermore at this point it's very critical to understand the fact that increasing $N$ alone and keeping $L$ fixed, does not provide a higher-resolution spectrum samples, it merely provides a smoother set of finely spaced samples of fixed DTFT.

The second relation is what gives the clues of expected spectral resolution of the DFT samples. Now from multiplication by a window, $w[n]$ of length $M$, property of the DTFT, it can be seen that a finite-length excerpt $v[n] = w[n] x[n]$ from a signal $x[n]$ of indefinite (or very long) length is providing you the DTFT spectrum of a limited-resolution version of the true spectrum $X(e^{j\omega})$ which is based on all the samples of $x[n]$

The limitation of the spectrum is fundamentally determined by the window length $M$ and window type. So in a practical setting, your FFT block length (or input buffer length etc) gives you the windowing effect. The effect can be roughly summarized as a blurring of the true spectrum. The result of which is that close-by (and detailed) frequency content is mixed up and get lost in each other, so that you cannot recover them back.

So if you want highly detailed spectral examination of a signal, you should observe it long enough (your window length should be long) and you should use a long enough DFT size $N$ so that you get a dense sample representation of the true (blurred) spectrum $X(e^{j\omega})$ after it was windowed in time domain whose effect is a convolution with a low-pass character frequency function $W(e^{j\omega})$ resulting in a smeared true spectrum.

Then coming to your practical setting: when you have a signal of $x[n]$ sampled at $F_s = 44.1$ kHz, and rectangularly windowed by $M=1024$ samples, the computed DFT samples are providing the values of blurred spectrum whose rought frequency reolution is about $2 F_s /1024 \approx 86.1$  Hz. Using a Hamming window would instead provide a resolution of 160 Hz. And you should expect a resolution something in between. Therefore, eventhough the DFT bins would sample the windowed DTFT at spacing of $f_k = F_s/N = 43$ Hz apart, you cannot expect meaningful information at less than $86$ Hz steps at best and $160$ Hz steps most usually.


displaying speech data in log scale using Matlab I have speech data that I can represent using matlab plot.
however, If i want to display it on log scale(like cooledit or Audition), I can take something like

data = audioread(pth);
for j=1:length(data)
    if data(j) >0
         data(j) = 20*log10(data(j));
    else
         data(j) = 20*log10(abs(data (j)));
    end
end

problem with this is that zero values map to infinity and negative values map to positive side . I need representation where Y axis would run in db like [0 -3 ..-21..-inf.. -21 -3 0](like in Cooledit) and the data is plotted as per this axis.
How do i do this?
sedy

It doesn't really make sense to display an instantaneous AC waveform on a log scale - as demonstrated by your results. I googled for cooledit screenshots and found this.

To accomplish this you would want to keep your data linearly scaled and draw it normally. Then just label the axis at the desired dB points.


what is the amplitude of quantization levels? okay , so there is an 8-bit linear PCM encoder . Now this mean that the amplitude of quantization levels is 2^8= 256 or I have to measure the values from oscilloscope ? does anybody have any idea how to answer that?

Practically speaking an 8-bit linear PCM encoder is implemented inside an 8-bit ADC (Analog to Digital Converter). Such an encoder accepts an analog input within a given refence voltage range and divides this range into 256 equally spaced quantization levels. Some minor deviations may occur based on the quantizer level distributions especially at the first and last levels.

At this stage analog input is converted into a PCM bit stream consisting of 8 bits of quantization codes per each sample.

Now to view these codes ? you can use a binary analyser to see each bit on or off in parallel...

Or if you mean observing the output of a DAC (Digital to Analog Converter) driven by this bit stream then you will be converting these bits back into analog and view it using the oscilloscope. The values will be steppy and be scaled by a reference voltage.


Why Is Super Resolution (SR) Possible? I've been reading about Super Resolution Image reconstruction (Reconstruction of high resolution image from multiple low resolution aliased images contain sub pixel shifts), and i want to know why SR is possible? this document explains why SR is possible, but i didn't understand completely!

please can you explain more clearly.

The more independent data you have, the more constrained are the possible solution sets that could produce that data, usually.  If any higher frequency content in the possible solution sets is constrained to not be completely arbitrary (which data derived from sub-pixel shifted sampling might so constrain), then the solution sets could possibly becomes sharper due to the steeper edges that can be reproduced using said higher frequency content.


Creating a time domain filter from a frequency domain plot In figure 11.23 of Monzingo's Introduction to Adaptive Arrays, the author uses a sinewave in the frequency domain to model the amplitude ripple in the passband of a channel of a receiver.

What does he mean by array bandwidth?
How does one create a time domain filter (FIR or IIR) with this custom/arbitrary frequency response?
What are other conventional channel amplitude and phase ripple filter models?

I don't know much about arrays, so I can only answer the second part of your question.

$$A(\omega)=(1+R\cos\omega T_0)\textrm{rect}\left(\frac{\omega}{2\pi B}\right)\tag{1}$$

where $\textrm{rect}()$ is the rectangular function. Equation (1) can be rewritten as

$$A(\omega)=\left(1+\frac{R}{2}e^{j\omega T_0}+\frac{R}{2}e^{-j\omega T_0}\right)\textrm{rect}\left(\frac{\omega}{2\pi B}\right)\tag{2}$$

which shows that $A(\omega)$ is a sum of three ideal low-pass filters with the same cut-off frequency, but with different delays ($0$, $-T_0$, and $T_0$) and different scaling factors. It is an ideal system and can only be approximated by a real system. The implementation suggested by (2) is a parallel connection of 3 low-pass filters with different delays and scaling factors. Note that in order for a causal filter to be able to approximate this ideal system reasonably well, you have to add a delay $\tau>T_0$. So the desired complex frequency response which you can base your design on is

$$H(\omega)=A(\omega)e^{-j\omega\tau}=\left(1+\frac{R}{2}e^{j\omega T_0}+\frac{R}{2}e^{-j\omega T_0}\right)e^{-j\omega\tau}\textrm{rect}\left(\frac{\omega}{2\pi B}\right)$$

Note that each of the 3 low-pass filters is a linear phase filter. So I would suggest you use FIR filters to approximate $H(\omega)$. You first design an FIR low-pass filter with the desired cut-off frequency with group delay $\tau-T_0$ (in samples). For a linear-phase filter, the delay determines the filter length $N=2*$(delay in samples)$+1$. So by choosing $\tau$ you also choose the filter length. If the impulse response of this filter is $h(n)$, and if $n_0$ is the number of samples corresponding to $T_0$, then the total impulse response is given by

$$g(n)=\frac{R}{2}h(n)+h(n-n_0)+\frac{R}{2}h(n-2n_0)$$


Notch filter: differences between IIR and FIR filters I'm trying to understand this great answer from Matt L. . It's said that "One advantage of IIR filters is that steep filters with high stopband attenuation can be realized with much fewer coefficients (and delays) than in the FIR case, i.e. they are computationally more efficient." First of all why this is true? Is this because of poles? Actually this comes from my previous question on removing 400Hz noise. There are two filters for that purpose. The first one is a FIR filter with following frequency response and pole-zero plot:

And the second one is a IIR filter:

Removing noise using the IIR filter was far better. Why this happens? According to the frequency responses, attenuation in the case of FIR filter is greater in the magnitude relative to IIR filter. Also I don't understand why zeros in the FIR filter are placed so that and how this positioning corresponds to that frequency response. I understand well the relation between pole-zero plot and frequency response for the IIR filter.

Your assumption why IIR filters can have steeper transitions from passbands to stopbands compared to FIR filters of the same order is correct: IIR filter have poles away from the origin of the complex plane, and poles inside the unit circle close to zeros on the unit circle cause the corresponding frequency response to change rapidly with frequency.
The FIR filter in your question has a much broader notch than the IIR filter. This is not only caused by the fundamental difference between FIR and IIR filters, but also because for some reason the FIR filter has two zeros close to the notch frequency instead of only one. This makes the filter more robust with respect to errors in the estimation of the noise frequency, but it may also attenuate desired frequency components. It's also much harder for an FIR filter to approximate a constant response away from the notch frequency because the frequency response of an FIR filter is a polynomial and not a rational function as is the case for IIR filters.
From the zero locations of the FIR filter you can see that the filter is a linear phase filter: the zeros are either on the unit circle or they are mirrored at the unit circle. Note that FIR filters don't need to have linear phase. The required filter order for a certain specification can often be somewhat reduced if we don't impose the linear phase constraint.


Evauating response to $x(n) = \sin{2\pi\frac{1}{4}n}$ given the system function Let the input signal $x(n) = \sin{2\pi\frac{1}{4}n}$ go in to the system described by: $$y(n) - y(n-1)+\frac{3}{16}y(n-2) =x(n).$$
  What is the output signal?

I've calculated the system function $H(z)$ to:
$$H(z) = \frac{1}{1-z^{-1}+\frac{3}{16}z^{-2}}.$$
But I don't know how to proceed. The book does the following:

$$H(w=2\pi\frac{1}{4}) = \frac{1}{1-e^{-i2\pi\frac{1}{4}} + \frac{3}{16}e^{-i2\frac{\pi}{4}2}} = 0.776e^{-i0.888}$$

I am not entirely sure where this comes from. Why can I put in the frequency in a function of $z$?

You can start from z-transform of sine:

$ Z[sin(2\pi f n)\cdot u(n)](z)= \frac{z \cdot sin(2 \pi f)}{z^2-2z \cdot cos(2 \pi f)+1}$

For your signal you have $f=1/4$ and substituting you obtain:

$ X(z)=\frac{z}{z^2+1} $

The output of the system will be:

$Y(Z)=X(Z) \cdot H(Z) = \frac{z}{z^2+1} \cdot \frac{1}{1-z^{-1}+\frac{3}{16}z^{-2}} $

Computing the inverse transform you will obtain:

$ y(n)=Z^{-1}[Y(z)](n)=\frac{2}{425} ((-64 + 52 i) (-i)^n - (64 + 52 i) i^n - 25 4^{(-n)} + 17 \cdot 3^{(n + 2)} 4^{(-n)})$

That I computed with wolfram alpha but you can solve it from yourself with partial fraction decomposition. Also if appears the $i$ coefficient the succession is not complex because is real for each integer $n \geq 0 $

The method used by your book uses the conversion from z-transform to dft ($ z=exp(j 2 \pi f)$ ):

$ X(z) \longrightarrow X(exp(j 2 \pi f)) $

The dft of digital sine is a train of Kronecker  Delta so the only part that remains doing the multiplication $ X(f) \cdot H(f) $ is $H(f)$ calcutated in $f=1/4+k$ and you can compute the idft of $H(exp(j 2 \pi (1/4+k))) $ to obtain the result.

But I think that its solution misses the $k$ that should be $\{0;1/2\}$ because the dft should be symmetric respect $f=0.5$ Nyquist frequency otherwise the idft begin complex. I don't know if the book defines an half dft for real signals. Also phase misses because sin has a phase of $\pi /2$ that it doesn't consider. I don't know exactly what is it doing, maybe it wants work with overlapping effects and the question misses other effects that need to be summed at this to perform the correct result.


Histogram evaluation for image quality concerning illumination I'm currently trying to develop an algorithm to control the camera characteristic of a camera (the camera characteristic can be controlled to adopt it to the different lightning situation which can occur)

In this case I need a good measurement of the exposure/lightning on the current image taken (i.e., if it's under/overexposed, well-exposed, maybe low/high contrast etc). And I'm therefore looking for a measurement on how to decide that. So far I couldn't find much about it. I found only ONE paper:

It's called "An optimal measure for camera focus and exposure" (DOI) from 2004 which describes some measurements like mean-values, entropy, standard deviation etc. based on the histogram to somehow measure the image quality (concerning focus and exposure). I also found many other papers which try to measure the image quality concerning things like blur, distortion etc etc but I don't need that obviously as I only want to control the camera characteristic to get a good contrast and exposure of my image -> thus, a well-balanced histogram.

Do you know any other measures which could help me with that or do you maybe have another idea on how solve that issue and how to build a good controller of my camera characteristic?

I understand that what you are looking for is a contrast measure. Since, I am not sure as to which different research papers you have looked at, I am giving here links to some important and good papers in this area. Refer to them. If you are unable to make sense of them, reply back, and I would try my best to help

Measuring pereptual contrast in color images

Measuring the perceived contrast of natural images

Hope this helps


How to achieve full number of words in a Bag-of-Features model based on Hierarchical K-Means? I'm implementing a baseline system that uses a Vocabulary Tree, i.e. a BoF model based on HKM, for images classification. Currently I'm obtaining low quality recognition due to the poor quality of the quantization structure resulting from the HKM.

I'm obtaining 100.000 words in my final vocabulary in a tree of depth 6 and branch factor 10 where the theoretical number of words is 10^6 = 1'000.000

In several papers, like those of Zisserman related to Large Scale Landmark recognition, they claim to be using a 1'000.000 words vocabulary, something I found difficult to understand since this number is theoretical while in practice there is no guarantee to obtain it.

Am I understanding something wrong? If not, what should I do to increase the vocabulary size despite of using more descriptors to train the tree?

PS: my only hint so far is using a different seeding algorithm for clusters initialization like k-means++ or gonzales

Thanks, any suggestion is welcome.

After reading @penelope's answer I came back to report my findings. The configuration of my implementation for building the vocabulary tree is: depth 6, branch factor 10, maximum 10 iterations, and random seeding.

Addressing some of the spots:

Severely underpopulated clusters at early levels: I checked clusters at first level and they were kind of uniformly distributed, e.g. the first cluster was typically  10% of the whole set.
K-means termination criteria: I'm running the clustering process until convergence or until reaching 10 iterations. I took this parameter from the experiments of Neister 2006 where the performance will not improve farther than ~12 iterations while doing so could affect the overall algorithm duration.
Uniform coverage of feature space: I am seeding centers using random selection, as I said in the question using other algorithm like k-means++ or gonzalez would become in a more uniform coverage of the feature space.
Number of training features should be considerably bigger than the expected number of clusters: with "expected number of cluster" I understand "expected number of words", being it so, this is very likely the root cause of my problem. I'm training a 1 million words vocabulary using 3 million descriptors while the experiments of Zisserman use ~17 million descriptors from the Oxford Buildings dataset and those of Marc Pollefeys team use ~16 million from the San Francisco Landmark dataset.
Pairing several detectors and descriptors: I have the restriction of comparing my approach which uses binary features against a baseline using DoG/MSER and SIFT descriptors. In the near future I plan to evaluate performance as well using other detectors like AGAST which produce more keypoints of similar quality to those produced by DoG.
HKM termination conditions: in @penelope's answer it was mentioned that obtaining a vocabulary tree of expected size was not a problem. I was wondering if such a thing was truly possible and started to think about HKM termination conditions. In my implementation they are:

Reaching last level.
Having less data than clusters.
Obtaining, at clustering initialization time, centers which are too close one another so there is a big chance they represent the same word, hence not very distinctive.

The first two conditions are enough and necessary but the third one isn't, I included it because the vocabulary might result more distinctive, at least in theory. I borrowed it from the hierarchical k-means index implementation in OpenCV. Could it be convenient removing it?

Unfortunately, I think what causes the problems is most probably the training set of images you're using to construct the tree. (or, more in general, the database you're working on). Without a way to inspect the concrete images, I can only give a few vague guesses.

First, I would diagnose where it goes wrong, by simply inspecting the tree and looking at the first branches which become severely underpopulated and consequently don't branch any more in the lower levels. I would look at those features and... think. But, since I can't actually do that since I don't have your dataset, let me try further.

If you're afraid that it's the K-means that is going wrong, not because the features aren't (more-or-less) uniformly distributed over the feature space, but simply because of the way that algorithm works, before implementing a fancier algorithm, you might first try running K-means several times for each branching and picking the best result (e.g. the one that divides the feature space in most equal parts). This will probably be too slow for the final implementation, but it could be a fast, minimum-effort way do diagnose the problem.

I would say, in order for HKM to work, your initial number of features should be several orders of magnitude greater than the expected number of clusters. So, if that's not the case, you know you're doing something wrong. I don't know how many is enough, but I would say minimum $100$ times, but better at least $1000$ times more features than expected clusters. This might be caused by two (related) issues:

Maybe your training set is just too small, or maybe it is a bad representation of the database, containing images that are too similar? Try increasing the size of the training set, and maybe (just for testing purposes), cherry pick the training set by selecting very different, representative images.
Maybe your feature detector performs badly on the image set you're using. For some images, some detectors just do not detect enough features. Problematic images usually contain large single-colored areas (e.g. white wall, blue sky). Try combining several feature detectors. (e.g. in my thesis I was working with DoG and MSER features). You can always do that, since you can finally describe all the features from several detectors with the same descriptor.

Unfortunately, there's no more ideas I can think of based on so little information. You might simply have an implementation error sometimes. I know it's not easy to here, and it's not a most helpful diagnosis, but as somebody who actually implemented such a CBIR system, I have to tell you, I never had any problem in constructing a vocabulary tree of expected size.

If you give any more info, or updates, or whatever, I can try and think of other ways this might have gone wrong. Until then, I hope this helps.


GnuRadio's SNR estimators usage I'm trying to calculate the SNR of a satellite signal using GRC, but all of the 4 offered estimators either have a lot of -NaN values or are ~35dB throughout the entire recording.
(I just connected the file source to the MPSK SNR estimator probe and redirected the messages to message debug.)

I've tried the examples adapted to my input and M2M4 estimator spits out somewhat expected values (so I've been told by others in the company).

The code:

from sys import argv
from scipy import fromfile, absolute, array_split, mean, sqrt, log10, ceil, linspace, complex64
from matplotlib import pyplot as plt

def snr_est_m2m4(signal):
    M2 = mean(absolute(signal)**2)
    M4 = mean(absolute(signal)**4)
    snr_rat = sqrt(2*M2*M2 - M4) / (M2 - sqrt(2*M2*M2 - M4))
    return 10.0*log10(snr_rat)

def chunkify(lst, chunk_size):
    return [lst[i * chunk_size : (i+1) * chunk_size] for i in range(int(ceil(len(lst)/chunk_size)))]

def main():
    if len(argv) != 2:
        return
    bits = fromfile(argv[1], dtype=complex64)
    m2m4_axis = [absolute(snr_est_m2m4(s)) for s in array_split(bits, 100000)]

    f1 = plt.figure(1)
    s1 = f1.add_subplot(1,1,1)

    chunk_size = 100
    smooth_m2m4 = [mean(c) for c in chunkify(m2m4_axis, chunk_size)]
    s1.plot(linspace(0, 1000, len(smooth_m2m4)), smooth_m2m4, 'b-', label='m2m4')
    s1.grid(True)
    s1.set_xlabel('time')
    s1.set_ylabel('snr')
    s1.legend()
    print("Done!")
    plt.show()

if __name__ == "__main__":
    main()

And the output I receive:

I know the recording is valid since the GRC chain we currently have implemented gives out an image of the Earth.
I've tried the estimators on a different signal but yet again, the values didn't make sense; the SNR was higher when there was no signal.

What do the parameters alpha and Samples between SNR messages in GRC represent? Why don't the estimators in GRC work as expected? What is the suggested way of calculating SNR?

I'm at a loss here, please send help.

The official documentation of that block:

  An SNR estimator for M-PSK signals that uses 2nd (M2) and 4th
  (M4) order moments. This estimator uses knowledge of the
  kurtosis of the signal ($k_a)$ and noise ($k_w$) to make its
  estimation. We use Beaulieu's approximations here to M-PSK
  signals and AWGN channels such that $k_a=1$ and $k_w=2$. These
  approximations significantly reduce the complexity of the
  calculations (and computations) required.

  Reference:
  D. R. Pauluzzi and N. C. Beaulieu, "A comparison of SNR
  estimation techniques for the AWGN channel," IEEE
  Trans. Communications, Vol. 48, No. 10, pp. 1681-1691, 2000.

The documentation of the constructor vaguely hints:

  alpha   the update rate of internal running average calculations."

I.e. probably a single-tap IIR ("exponentially weighted moving average") with feeback coefficient alpha.

Verifying that with the source code:

for(int i = 0; i < noutput_items; i++) {
    double y1 = abs(input[i])*abs(input[i]);
    d_y1 = d_alpha*y1 + d_beta*d_y1;

    double y2 = abs(input[i])*abs(input[i])*abs(input[i])*abs(input[i]);
    d_y2 = d_alpha*y2 + d_beta*d_y2;
}

Yes, that alpha is just the coefficient of a single-tap IIR (further up you'll find d_beta = 1 - d_alpha).

So, this estimator should do exactly what it promises: for a zero-mean signal (where the variance is simply $\text{Var }X = \text E(X^2)$), it calculates the second and fourth central moment and then just

double
mpsk_snr_est_m2m4::snr()
{
  double y1_2 = d_y1*d_y1;
  d_signal = sqrt(2*y1_2 - d_y2);
  d_noise = d_y1 - sqrt(2*y1_2 - d_y2);
  return 10.0*log10(d_signal / d_noise);
}

So, it does the usual kurtosis-based estimate.

By the way, I'd fully agree the parameter documentation – to put it friendly – leaves a bit of room for improvement. Hit us / me up if you'd be willing to contribute a bit of a documentation improvement!


Where is the mistake in the convolution? Let :
$$
x[n]=\begin{cases}
1&\text{if $0\leq n\leq 4$}\\ \\
0&\text{if otherwise}
\end{cases}
\qquad
\text{and}
\qquad
h[n]=\begin{cases}
\alpha^{n}&\text{if $0\leq n\leq 6$}\\ \\
0&\text{if otherwise}
\end{cases}
$$
We must compute the convolution :
$$
(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]
$$
Observe that $x[n]=u[n]$ for $n\in[0,4]$ and $0$ otherwise. Furthermore, $h[n]=\alpha^{n}u[n]$ for $n\in[0,6]$. Thus :
\begin{array}{c | c c c c c c c }
\mathbf{x} & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
\hline
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
0 & 0 & 0 & 0& 0& 0 & 0 & 0 \\
0 & 0 & 0 & 0& 0& 0 & 0 & 0
\end{array}
\begin{align*}
\displaystyle (x*h)[n]&=\left[\sum_{k=0}^{0}\alpha^{k}\;,\;\sum_{k=0}^{1}\alpha^{k}\;,\;\cdots\;,\;\sum_{k=0}^{4}\alpha^{k}\;,\;\sum_{k=1}^{5}\alpha^{k}\;,\;\sum_{k=2}^{6}\alpha^{k}\;,\;\cdots\;,\;\sum_{k=5}^{6}\alpha^{k}\;,\;\sum_{k=6}^{6}\alpha^{k}\right]\\ \\
&=\begin{cases}
\displaystyle\sum_{k=0}^{n}\alpha^{k}&\text{if $0\leq n\leq 4$}\\
\displaystyle\sum_{k=1}^{n}\alpha^{k}&\text{if $n=5$}\\
\displaystyle\sum_{k=n+2}^{6}\alpha^{k}&\text{if $0\leq n\leq4$}
\end{cases}
\end{align*}
I wish to know where my mistake is and how to fix it. Thank you :)

$3^{rd}$ term in your final answer should have $6 \le n \le 10$, instead of $0 \le n \le 4$. And, the limit of summation will go from $n-4$ to $6$.
$$\sum^{6}_{k = n-4} \alpha^k \ ,  6 \le n \le 10$$
$0 \le n \le 4$ case has already been covered by the $1^{st}$ term of final answer.


fftshift in MATLAB with even number of data points in double sided spectrum I have a question with reference to this Table.

With even N, the frequency axis extremes should be $\pm$Fs/2, where Fs is the sampling frequency. However in the array we have only one value corresponding to the Nyquist frequency. Following the notation of the table in the first column, N=10, we have

[$c_{0}$ $c_{1}$ $c_{2}$ $c_{3}$ $c_{4}$ $c_{-5}$ $c_{-4}$ $c_{-3}$ $c_{-2}$ $c_{-1}$]

which after fftshift should become

[$c_{-5}$ $c_{-4}$ $c_{-3}$ $c_{-2}$ $c_{-1}$ $c_{0}$ $c_{1}$ $c_{2}$ $c_{3}$ $c_{4}$  ]

The $c_{-5}$  value corresponds to the Nyquist frequency. So how is our double sided frequency spectrum symmetric when N is even? For real valued function it will be zero, but is there a reason that MATLAB calculates the negative frequency first as shown in the Table?

Example:
A_even=[ 0 0 0 1 1 1 1  0 0 0],

N=10

B=fft(A_even)'

C=fftshift(B)

Then C is equal to

$$\begin{matrix} 0.0000 + 0.0000i\\
-0.3090 - 0.9511i\\
0.4271 + 0.5878i\\
0.8090 + 0.5878i\\
-2.9271 - 0.9511i\\
4.0000 + 0.0000i\\
-2.9271 + 0.9511i\\
0.8090 - 0.5878i\\
0.4271 - 0.5878i\\
-0.3090 + 0.9511i\\
\end{matrix}$$

This implies $c_{-5}$ =  0.0000 + 0.0000i  and $c_{0}$= 4.0000+0000i
Thanks.

If the size $N$ of the DFT is even, only one "extremal" point (after fftshift) is Nyquist.  If $N$ is even, you cannot have an arbitary $c_{-5}$ and $c_5$.  They must add to be whatever your $c_{-5}$ term is.

If the input to the DFT is purely real and if $N$ is even, consider the Nyquist point, at $c_{N/2}$ to be split in half.  One half is the negative-frequency component and the other half is the positive-frequency component.  Since the negative and positive frequency components must be complex conjugates of each other, and since $c_{-N/2}=c_{N/2}$ after splitting, that means that $c_{N/2}$ must be purely real just as $c_0$ must also be real.


Confusion in initial condition of differential equation using Laplace transform transform I'm confused in solving linear constant coefficients differential equations (LCCDEs) by Laplace transform if initial conditions are given at time

just before $t=0$
just after $t=0$
exactly at $t=0$

Is method of solving LCCDE by Laplace transform is same in all three cases or it's different?
I know it's more of mathematical question and I asked this question in mathematics stack exchange but I didn't get answers, so I ask it here as it is equally a problem of systems differential equation.!

Initial conditions are always given at $t=0^-$, because they define the state of the system before any input is applied, and - by definition - the input is applied at $t=0$. The state at $t=0^+$ is determined by the initial conditions as well as by the input signal.
The unilateral Laplace transform can be used to solve LCCDEs with initial conditions $y(0^-), y'(0^-),\ldots$ because of the definition
$$\mathcal{L}\{f(t)\}=F(s)=\int_{0^{\color{red} -}}^{\infty}f(t)e^{-st}dt\tag{1}$$
from which it follows that
$$\mathcal{L}\{f'(t)\}=sF(s)-y(0^-)\tag{2}$$
Note that it's common to write initial conditions as $y(0),y'(0),\ldots$, when actually $t=0^{-}$ is meant.
EDIT: Concerning the example in the book: if there is no Dirac delta impulse in the current $i(t)$ at $t=0$, the capacitor voltage $v_c(t)$ cannot jump at $t=0$. Consequently, $v_c(0^-)=v_c(0)=v_c(0^+)$ must hold.
In general, the initial conditions define values of the output signal and its derivatives right before the source signal is switched on. It can be the case that the limits of the output and its derivatives exist at $t=0$. If that is the case, it doesn't make a difference if we use $t=0^{-}$ or $t=0$ or $t=0^{+}$ because the corresponding function values are all the same. If the output signal or its derivatives are discontinuous at $t=0$, then the distinction becomes important, and the initial conditions define the values just before the discontinuity, i.e., at $t=0^{-}$.


creating a seamless signal / loop using interpolation I'm trying to create a seamless loop using a "non-periodic" signal using interpolation to smooth out the beginning and the end but I'm still getting a click at the beginning when it loops and I listen to it.  Can interpolation do this?

Please Note that the equation is done this way to create a "non-periodic" signal to test when a signal does not start and end at the zero crossings.

I used audacity to zoom in and play the loops.  I've included images along with the matlab/octave code.

% combines sig to create seamless loop
clear all,clc,tic;
fs=44100;
dirpathtmp=strcat('/tmp/'); %/home/rat/Documents/octave/eq_research/main/transform/voice
fn=strcat('atest1',sprintf('%02d',fs),'.wav');

t1=linspace(0,2*pi,fs); %need to round off no decimals
t2=linspace(0,2*pi,fs); %need to round off no decimals
t3=linspace(0,2*pi,fs); %need to round off no decimals

%Create signal in different arrays
y1=sin(7*t1);
y2=sin(12*t2);
y3=cos(2.2*t3);
yt=y1+y2+y3;

ytnorm=(yt/max(abs(yt))*.8); %normalize signal

%change end points to stddev number
ytnorm_std=std(ytnorm);

%for loop to replace points to create seamless loop
ytnorm(1)=0; %set first point to 0
ytnorm(fs)=0; %set last point to zero

for kk=2:1:15;
    kk
    ytnorm(kk)=interp1(t1,ytnorm,kk,'spline');
    ytnorm((fs+1)-kk)=interp1(t1,ytnorm,((fs+1)-kk),'spline');
end;

wavwrite([ytnorm'] ,fs,16,strcat(dirpathtmp,fn));

plot(ytnorm)

fprintf('\nfinally Done-elapsed time -%4.4fsec- or -%4.4fmins- or -%4.4fhours-\n',toc,toc/60,toc/3600);

An alternate method is to create an array of two copies of your input signal

y = [ x , x ]         % MATLAB
y = hstack(( x , x )) # NUMPY

Run $x[k]$ through a LPF.  Then slice $N$ samples from the "middle" of $y[k]$ where $N$ is the length of $x[k]$.

xsmooth = y( o : o + N - 1) % MATLAB
xsmooth = y[ o : o + N ]    # NUMPY

Where $o$ is an appropriately chosen offset.  Be mindful of filter "tails" if you carry out the LPF through conventional convolution.


How to interpret the window size and FFT size while performing FFT? I have a discrete signal of length 98,000 samples and I am supposed to a FFT of the discrete signal with certain window size and FFT size. Can someone help me explain what is meant by the window and the FFT size here?
Also, what would the size of the output of the signal post FFT?

Thanks in advance!

There are lots and lots  of related questions on this topic in this SE and also on many discrete signal processing books.

https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf
https://holometer.fnal.gov/GH_FFT.pdf


Fringe movement analysis to extract displacement I am using a Michelson interferometer with the interferogram recorded on a CMOS sensor. I have taken a video of the fringes moving when a displacement of 50 $\mu m$ is imposed to one of the mirrors.
The frame rate is 17.9980 and it is already the maximum I can have.
For a 50 $\mu m$ displacement I should see around 79 fringes passing on the sensor, as each fringe that passes corresponds to one full wavelength (632.8 nm). However, I cannot obtain clear "peaks" that allow me to identify if a fringe is moving.
In the picture you can see the frame with the fringes (top left) and the signal in the central pixel column (top right) for each frame. Then, I transformed the signals using FFT in Matlab, and I obtained the spatial frequency (always the same, so it should be correct) and the phase. Theoretically, I should be able to extract the displacement in two consecutive frames checking the phase difference, so that $$\Delta x = \frac{\lambda}{2} \frac{1}{2\pi}\Delta\phi$$
and the total displacement being the sum of the $\Delta x$ for each couple of frames.
Alternatively, I thought that I could extract the displacements looking at the integral of the signal $$\Delta S_k = S_k - S_0 $$, which is the difference between the signal at the $k$-frame and the initial signal at frame $0$. When the two signals are in phase, the integral of the difference signal should be close to zero, while when the two signals are not in phase the integral increase to a maximum when the $\Delta \phi=\pi$ (bottom right)

In the following picture there is the wrapped phase (top), unwrapped phase (centre) at the main frequency, and the integral of the difference signal at each frame (bottom). Imposing a 50 $\mu m$ in the same direction, the sign of $Delta \phi$ should be always the same, but apparently it is not. As you can see, from the difference signal the peaks are not easy to determine, there are "double peaks" and small peaks. Theoretically, I should find 79 peaks but I got more.

Is there an easier way to count fringes?
I guess that the peaks would be easier to identify if I move the mirrors more slowly?

Answer to: Is there an easier way to count fringes?
I synthesized an image to resemble your image with
x = linspace(0, 1, 256);
[x,y] = meshgrid(x,x);
image = sin((7*x + 2*y) * 2*pi) ./ (0.1 + (x - 0.5).^2 + (y - 0.6).^2);
imagesc(image);

Then I can get the vector the 2D FFT
A = abs(fft2(image));
A = A(1:end/2, 1:end/2)
[~, imax] = max(abs(A(:))) % 899
% get line and column from the index
idx = [floor((imax-1)/size(A, 2))+1, mod(imax-1, size(A, 2))]

idx is set to the expected wave vector [7, 2] thus the length of the fringes must be $l = \sqrt{W/7^2 + H/2^2}$, where W and H are width and height of the image.
Update after comment

What kind of information can the wave vector provide about the fringe movement? In my video the interference pattern does not change, it just moves...

You must have realized that with this interference pattern you can only detect the component in the direction of the wave vector.
Also you cannot detect movements beyond one half of the distance between the fringes, so you can use unwrap to solve it. Here a demonstration with multiple frames and some movement.
function fringes
  k = [3, 10];
  t = linspace(0, 1, 100);
  phi = t * 0;
  d = 2*sin(5*t);
  WATCH = 0; % set to 1 if you want to see all frames
  % simulate some movement on the fringes
  for i = 1:length(t)
    image = create_image(k, d(i), 256);
    %% if you want to watch the movement
    if(i == 1 | WATCH)
      figure(1), imagesc(image), drawnow;
    end
    [~, phi(i)] = detect_fringes(image);
  end

  figure(2);

  dist = sqrt(256/k(1)^2 + 256/k(2)^2);
  c = plot(t, d*dist, ...
       t, phi*dist/pi/2, ...
       t, unwrap(phi)*dist/pi/2 ...
  );
  l = legend('displacement', 'pricipal phase', 'estimated displacement')
  set(l, 'fontsize', 12);
  set(c, 'linewidth', 2);
end

function image = create_image(k, d, N);
  x = linspace(0, 1, N);
  [x,y] = meshgrid(x,x);
  image = sin((k(1)*x + k(2)*y + d) * 2*pi) ./ (0.1 + (x - 0.5).^2 + (y - 0.6).^2);
  image = image + randn(size(image)) * 10;
end

function [k, phi] = detect_fringes(image)
  A = fft2(image);
  A = A(1:end/2, 1:end/2);
  [~, imax] = max(abs(A(:)));
  % get line and column from the index
  k = [floor((imax-1)/size(A, 2)), mod(imax-1, size(A, 2))];
  phi = angle(A(imax));
end

I added some noise to make it more realistic

The result is this

Tested in Octave 6.1.0
Edit: Why the wave length is inversely propoertional to |k|
The wave length is the separation between two points where the wave field is in phase.
The wave vector $\vec{k}$ is associated with a wave of the form
$$
\begin{eqnarray}
 F &=& A exp\left( 2i \pi \vec{k} \cdot \vec{u}\right) \\
   &=& A exp\left( 2i \pi |\vec{k}| \cdot |\vec{u}| cos( \angle(\vec{k}, \vec{u}))\right) \\
\end{eqnarray}
$$
Let $\phi$ be the phase of the field
$$ \phi(u) = |\vec{k}| \cdot |\vec{u}| cos( \angle(\vec{k}, \vec{u})) $$
One wave length may be defined by the minimum $\Delta \vec u$ such that $\phi(u) + 2 \pi = \phi(u + \Delta \vec u)$, the direction of $\Delta \vec u$ is choosen to maximize $cos(\angle(\vec k, \Delta \vec u))$, i.e. $\angle(\vec k, \Delta \vec u) = 0$, thus $\Delta \vec u$ is parallel to $\vec k$. In this situation $\phi(u + \Delta \vec u) = \phi(u) + |\vec k| |\Delta \vec u| = \phi(u) + 2\pi$. Solving for $\lambda = |\Delta \vec u| = 2\pi / |\vec k|$.


FFT filter to remove cd scratch I have a .wav file which came off of a cd with a scratch, so there is a periodic "clicking" of about 6 Hz.  I was playing around with trying to squash the clicking by taking FFT and filtering the range around this frequency, but with no success.  I tried smoothing the low end of the FFT spectrum, but it didn't change anything.  Then I tried zeroing this part of the spectrum up to higher and higher indices, but it still wouldn't kill the click.  Eventually it started to attenuate the entire signal, but the click was still there, and perhaps even more prominent than before (relative to the rest of the signal).  Why is this the case?

My first speculation is that the click is somehow not quite steady enough to be localized in a single region of the spectrum.  Can anyone confirm this hypothesis or else explain the true source of the issue?  Also, is there a better means of eliminating the click?

This approach with notch filter with not work. All clicks are impulse-like sounds and we know that an impulse has frequency content at almost every frequency. What you are trying to do, by applying the $6 \ \texttt{Hz}$ notch filter is to remove sinusoid with that particular frequency in your signal, isn't it? So do you hear that your clicks are such low-frequency sounds? Definitely not! Those clicks are spreading into whole of your spectrum (like very short burst of noise), and filtering out $6 \ \texttt{Hz}$ won't help!.
You also didn't shared any more informations about nature of your clicks and a signal (spectrogram, waveform or even recording would be great), because it might be the case that your signal is also impulse-like (i.e. recording of a shooting range), which will make the recognition and cleaning task bit more complicated. Also are your clicks the same over time? This could be also the case for digital systems, and simple correlation calculation should do the job. Additionally if you can sacrifice some of your recording quality and click's are very short, then you can use the median filtering. You can also use some time-frequency methods (based on wavelets, spectrogram calculation, LPC's) or  statistical approach, such as AR modelling, etc.
Here is some reading for you - I hope it will guide you and you will find appropriate method for your kind of signals.

Audio Click Removal using LPC
Godsil S., Rayner P. - Digital Audio Restoration - lot's of references
AR modelling based approach - book chapter
Detection of clicks using sinusoidal modelling for the confirmation of the clicks
Automatic detection and removal of impulsive noise in audio signals


FIR filter design with matlab firls() function I have defined the following linear phase filter of order 4 (real and symmetric) in matlab:
h_d = [0.8367 -1.537 1 -1.537 0.8367]
I calculate the frequency response (for many points) and then I use firls function (which is supposed to give the least-squares solution) to find a linear phase FIR filter (5 taps) that matches the desired response. However the filter is nowhere close to h_d. What am I doing wrong?
Here is the code:
f = [0:1/(2^14-1):1];
[a,~] = freqz([0.8367 -1.537 1 -1.537 0.8367],1,f*pi);
a=abs(a);
n = 4; % Filter order
b_lp = firls(n,f,a);
[h,w] = freqz(b_lp,1,512);
figure(1)
plot(f,20*log10(abs(a)),'-s')
hold on
plot(w/pi,20*log10(abs(h)))

I appreciate any insights.

The function firls() is meant to design filters with piecewise constant magnitude responses. So in practice you use only a few frequency points and the corresponding desired magnitude values, and the function computes a linear interpolation between the given frequency points. Of course, in theory your call to firls is correct, but I guess that the resulting system of linear equations becomes ill-conditioned.
I wrote a function lslevin.m, which can be used in the way you intended to use firls:

h = [0.8367 -1.537 1 -1.537 0.8367];
[H,w] = freqz(h,1,2048);    % that's more than enough frequency points
h2 = lslevin(5,w,H,ones(size(H)));
[H2,w] = freqz(h2,1,2048);
plot(w/pi,abs(H),w/pi,abs(H2),'r--')


What is the meaning of the DFT? Possible Duplicate:
  Real Discrete Fourier Transform
  What is the most lucid, intuitive explanation for the various FTs - CFT, DFT, DTFT and the Fourier Series?
  Discrete-time Fourier transform

I read about the Discrete Fourier Transform. I understand how to compute it given a sequence of numbers(digital signal), but I still didn't get what's the meaning of the new sequence I get as the output.

I mean, I know that the Fourier Transform takes an analog signal representation in the time domain and returns its representation in the frequency domain, its spectrum.
If I look at the analog signal in a time-window $T$, sampling it $N$ times , then the sampling interval will be $t_0 = T/N$ and the sampling frequency will be $f_s =1/t_0$.
The sampled values of the analog signal give me a new digital signal, lets say $s_n$

So what is the meaning of the DFT of this digital signal $s_n$? How the DFT sequence related to $s_n$? How it is related to the original analog signal? What is the meaning of these frequency components I get? Should they be taught of an approximation the previous analog signal frequency components or as sampled values of its spectrum?

Thanks a lot!

As you noted, the discrete Fourier transform (DFT) maps a length-$N$ sequence to its (also length-$N$) frequency-domain equivalent. The time-domain signal $s[n]$ is transformed into a frequency-domain signal $S[k]$, such that the following two relations hold true:

$$
S[k] = \sum_{n=0}^{N-1} s[n] e^{\frac{-j2\pi k n}{N}} \text{ (the DFT)}
$$
$$
s[n] = \sum_{k=0}^{N-1} S[k] e^{\frac{j2\pi k n}{N}} \text{ (the inverse DFT)}
$$
(ignoring scale factors that can be placed in many different places depending upon the preference of the author; one way to make them "symmetric" is to place a factor of $\frac{1}{\sqrt{N}}$ in front of each equation: the unitary DFT).

There are a number of intuitive interpretations for the action implied by these equations. Here's the way I usually think about them:

Forward DFT: Each DFT output value $S[k]$ is a complex value that represents the amplitude and phase of $s[n]$'s content at frequency $\frac{2\pi k}{N}$. The multiplication by the exponential function effectively shifts $s[n]$ down in frequency by $\frac{2\pi k}{N}$. The sum over $N$ time samples can then be thought of as applying a decimating lowpass filter. So, effectively, the $N$ outputs of the DFT represent the results of applying a bank of equally-spaced filters across $s[n]$'s frequency band, thus measuring "how much" energy is present in various frequency bands in the input signal.
Inverse DFT: The DFT output values $S[k]$ are used to weight a number of complex exponential functions in order to resynthesize the original time-domain signal $s[n]$. The frequency associated with $S[k]$ is the additive inverse of the frequency used during the forward DFT. This view clearly shows that the DFT can be thought of as a means to decompose an arbitary finite-length time-domain signal $s[n]$ into a set of weighted orthogonal complex exponential functions (or "complex sinusoids").

Now, extending this concept a bit: when you sample a continuous-time signal $s(t)$ at discrete time instants, you get a (potentially infinite-length) discrete-time signal $s[n]$. Applying the Fourier transform to such a discrete-time signal results in the discrete time Fourier transform (DTFT), which is continuous in frequency and, like the DFT, periodic in frequency with period $2\pi$. Assuming that a large-enough sample rate was used during the conversion to discrete-time, the result of the DTFT will look a lot like the Fourier transform of the original signal $s(t)$.The DFT is only applicable to finite-length signals, so in order to apply it, one must truncate the discrete-time signal $s[n]$ to some finite length $N$, resulting in a potentially-shorter signal $s_N[n]$.

For this special case of a finite-length discrete-time signal, the $N$ DFT outputs are just equally-spaced samples of $s_N[n]$'s DTFT.
$s_N[n]$'s DTFT is related to $s[n]$'s DTFT; they may differ if $s[n]$ was longer than $N$ samples long originally and some truncation/windowing was involved in generating a finite sample length.
$s[n]$'s DTFT is related to the original signal's Fourier transform, with the caveat that the DTFT is periodic, so all of the frequency content from the continuous-time signal that is unambiguously representable in discrete-time (the bandwidth of which is determined by the sample rate used) is contained within a single length-$2\pi$ period of $s[n]$'s DTFT.


Does audio upsampling create noise/artifacts or degrade the signal? I've read that upsampling performed in digital music playback can color the sound, produce artifacts, etc.  For example, an audio file ripped from a CD might be 44.1Khz/16-bit, and then upconverted to 48Khz/16-bit and played via an optical digital audio output.  Audiophiles say this is bad because the signal needs to be "bit perfect" to be reproduced correctly.

Is this correct?  My vague knowledge of DSP from grad school leads me to think that all upsampling should do is increase the bandwidth of the signal.  But, since the source signal is bandlimited, there shouldn't be anything new in the added bandwidth.  And I don't see why the process would color the sound.

What am I missing?

The up-sampling process will always change the signal in some measurable way. However, if it's done properly the changes are negligible and don't result it any audible difference. Most commercially sample rate converters (hardware or software implementations), do a really good job at this.

Off course, if done badly, upsampling can result in clearly audible signal degradation. I'm not familiar with Apple's implementation but I would assume that they got this correctly.


Generate the Convolution Matrix of 2D Kernel for Convolution Shape of `same` I want to find a convolution matrix for a certain 2D kernel $ H $.
For example, for image Img of size $ m \times n $ , I want (in MATALB):

T * Img = reshape(conv2(Img, H, 'same'), [], 1);

Where T is the convolution matrix and same means the Convolution Shape (Output Size) matched the input size.

Theoretically, H should be converted to a toeplitz matrix, I'm using the MATLAB function convmtx2():

T = convmtx2(H, m, n);

Yet T is of size $ (m+2) (n+2) \times (mn) $ as MATLAB's convmtx2 generates a convolution matrix which matches Convolution Shape of full.

Is there a way to generate the Convolution Matrix which matches using conv2() with the same convolution shape parameter?

I cannot test this on my computer because I do not have the convtmx2 function, here is what the MATLAB help says:

http://www.mathworks.com/help/toolbox/images/ref/convmtx2.html

T = convmtx2(H,m,n) returns the convolution matrix T for the matrix H. If X is an m-by-n matrix, then reshape(T*X(:),size(H)+[m n]-1) is the same as conv2(X,H).

This would get the same resulting convolution of conv2(X,H) but then you would still have to pull out the correct piece of the convolution.