Linear approximation of a system described by a logarithm - how? The system described by a logarithm of any base, let's call it $y(t) = \ln(t)$, is non-linear.

Is it possible, by any means (transforms, operators etc.), and for small values of $t$ (for example, $t>0$ and $t\leq 0.5$), to obtain a linear approximation of the system (in such $y(t_1+t_2) = y(t_1)+y(t_2)$) with an acceptable percentage error (i.e, up to $5\%$)?
Is there any log base more suited to perform the task?

Functions like $\sqrt{\cdot}$ or $\log{(1+\cdot)}$ may emulate nonlinear system behavior like saturation. They are although used as companding methods, such as the $\mu$-law (or the A-law):
$$ x\to \mathrm{sgn}x\frac{\log{(1+\mu|x|})}{\log{(1+\mu)}}\,. $$
This can be used prior to quantization for signals with high-dynamics, and allows one to use uniform quantization and maintain a similar amount of relative errors for weak as well as strong signals (see for instance Digital filtering using logarithmic arithmetic (I glad to discover this work by Nick Kingsbury).
For telephony, people developed segment ("piecewise linear", but actually piecewise affine to be correct) companding versions of such laws, for instance H. Kaneko, "A Unified Formulation of Segment Companding Laws and Synthesis of Codecs and Digital Compandors". Basically, $\log(x+y)$ can be approximated by some affine law in its proper interval. I am aware of such approximations  for square roots: if $a\ge b$, then for instance with about 4% precision:
$$\sqrt{a^2+b^2} \sim 0.96a+0.4b  $$
yet I have to check for logarithms. In the same framework, as already said by @ben, if $a$ is a fixed quantity:

if $a>b$, then $\log(a+b)\sim b/a+\log a$

which is  an affine approximation.


Demodulation of 4-QAM Any help to get me started is appreciated. I have an array of values for the 4-QAM where attach is a copy of part of the values of the array. I want to know how do I get the bits 0 and 1s from those values. I would be using python to program. Thank you very much

I have searched the forum and the post below is something that I am trying to do..I understand on how to get the bits. Is there a formula where I can feed the above values?  Most references state to check the BER (bit error rate) but do we get the bits..

Recover IQ samples from signal

"If you just want to demodulate, then it is not necessary to recover I and Q components. What you have is most likely the QAM symbols. For example, for 16-QAM they are 16 different integer values that are assigned to different constellation points. You need to know the type of assignment (i.e. coding) that is used in the QAM constellation at the modulator."

It is not completely clear what is the question but based on the last quoted paragraph, if you need to know how the 4-QAM mapping was done at transmitter. For example, if the mapping is

bits {00} -> 0

bits {01} -> 1

bits {10} -> 2

bits {11} -> 3

then if you receive an integer stream $\{2,3,0,1,2,1\}$, then the bits are $\{10,11,00,01,10,01\}$. Does that answer? If not, what exactly is your input array?


Choosing clock frequency and number of bits of an OFDM transceiver I want to implement an OFDM transceiver chip and my question is how to choose the clock frequency of the IFFT block? The number of output bits of the IFFT block? Also, how has the clock frequency and the number of output bits of the digital-to-analog converter (DAC) to be chosen?

The clock frequency $f_\mathrm{s}$ of the digital-to-analog converter (DAC) ultimately depends on the desired bitrate $R$ and some other parameters like guard interval length and number of bits per subcarrier. You can find the detailed formulas in a related question.

The speed of the IFFT block depends on $f_\mathrm{s}$ and the number of subcarriers $N_\mathrm{SC}$. It has to transfrom $f_\mathrm{S}/N_\mathrm{SC}$ sets of values each set consisting of $N_\mathrm{SC}$ complex numbers. What clock frequencies are involved in such an IFFT block depends on implementation details (mainly degree of parallelization) and cannot be answered in general.


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

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

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

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

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

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

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

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

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

Do you have any suggestions?

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

change

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

in your codes to

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

You should be able to get the right skeleton.


DFT's Unitary Property I have doubt, please clear.

A transform is said to be unitary transform when the product of the transform matrix and its conjugate transpose (Hermitian) matrix is equal to Identity matrix.

DFT for N= 4 gives

c =

 4     0     0     0
 0     4     0     0
 0     0     4     0
 0     0     0     4

Thus, it satisfies the Unitary property

But when N=2, we have

2      1+j
1-j     2

This doesn't satisfy unitary property right? Then how can we say the DFT is a Unitary Transform? Where I have misunderstood the concept? Help

Just so we're on the same page: given some sequence $x[n]$ for $n=0, \dots, N-1$, the $k$'th DFT coefficient is (forgetting the scaling factor)

$$X[k] = \sum_{n=0}^{N-1} x[n] \exp(-j 2\pi k n/N)$$

We can get all the coefficients by writing this as $ X = Dx$, where $D$ is the $N \times N$ DFT matrix and $D_{kn} = e^{-j 2\pi k n /N}$. For the N = 2 case,

$$D = \begin{pmatrix}
e^{-j 2\pi (0)(0) /2} &  e^{-j 2\pi (0)(1)/2}\\
e^{-j 2\pi (1)(0) /2} &  e^{-j 2\pi (1)(1)/2}
\end{pmatrix} =
\begin{pmatrix}
e^{-j 0} &  e^{-j 0}\\
e^{-j 0} &  e^{-j \pi}
\end{pmatrix} = \begin{pmatrix}
1 &  1\\
1 &  -1
\end{pmatrix},$$
for which $D^* D = D^T D = D^2 = 2I$, so (minus the scaling factor), unitary.


What are the audio equivalents of images like "Lena", "Mandrill", and "Cameraman"? I don't know what to call these images, but they're very often used in the image processing world as standard examples to show the output of a particular operation. Here's a link to the page describing Lena as a "standard test image". Are there similar canonical audio clips used to demonstrate the effects of certain audio filters or operations?

The closest example I can think of is the beginning of Suzanne Vega's "Tom's Diner" which has been used for the mpeg-1 layer 3 development, and is still occasionally used to demo audio codecs.


Why edge sharpening produces high frequency? I have a low-resolution image in which the high frequencies are missing. When I apply an edge sharpening filter some of the missed high frequency is recovered.

I am wondering why this edge sharpening process in the image domain produces new high-frequency information in the frequency domain? Is this new information reliable?

Edges are not the best defined features in images.

However, they can be associated, locally, at a certain scale, with relative variations in intensity along a first direction, combined with a relative smoothness in a complementary (for instance orthogonal) second direction. When looking along the first direction, the 1D intensity profile exhibit variations, faster  than in the other direction. Several typical edges profiles exist, as depicted below.

Those are typically enhanced by derivative-like operators, which are weakly sensitive to flat, linear or generally regular intensity variations. In other words, they behave as band-pass or high-pass operators.
 In the other direction, gradient detectors offer more smoothing (more low-pass).

Another vision is to smooth an image, and subtract this result  from the original image. By removing low-frequencies, higher ones are relatively amplified. But to avoid noise amplification, some smoothing is often included in the process.

Hence, linear or not, edge operators tend to emphasize mid- to high frequencies relatively to the lower ones. However, all this may happen at different scales, depending on the edgy objects of the picture.

Now,  Is this new information reliable?

Not really, if you rely on mere filtering or enhancement. It is easy to produce spurious edges. But within a restoration or inverse problem context, with models and an objective function to minimize, more reliable information can be obtained, as long as the models mimic well the physics of your acquisition scheme.


Appying a Filter Several Times on Data I have data from an experiment and apply a lowpass filter (digital Butterworth filter, in Matlab).

My question ist what happens if I apply the same filter e. g. 2 times on the data?

What is changing? Is it like I have the double order?

The question refers to the filtfilt function in Matlab:

http://www.mathworks.de/de/help/signal/ref/filtfilt.html

The source [1] says on page 336:

  In many filtering problems, we would prefer that the phase
  characteristics be zero or linear. For causal filters, it is
  impossible to have zero phase. However, for many filtering
  applications, it is not necessary that the impulse response of the
  filter be zero for n < 0 if the processing is not to be carried out in
  real time. One technique commonly used in discrete-time filtering when
  the data to be filtered are of finite duration and are stored, for
  example, in computer memory is to process the data forward and then
  backward through the same filter.

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

How about answer your question in an ultimate way? What will happens to if you apply a low pass filter infinite times or practically a large number of times? I believe it is more meaningful for you to understand this one.

I want to first emphasize the fact that apply a filter $F$ to $d$ for $k$ times, is equivalent to apply a filter $F'$ to $d$ for once, if we have
$$F' = \underbrace{F*F*\cdots*F}_{k \,{\rm times}} $$.

The answer to the 'ultimate' version of your question is that: no matter where you start with your lowpass filter, say $F=[f_0,f_1,\cdots,f_n]$, with $\sum_{i=0}^{n} f_i = 1$ and $f_i \geq 0$, you will end up to apply a Gaussian filter $G$ when you apply filter $F$ to some discrete signal $d=[d_0,d_1,\cdots,d_m]$ for $k\rightarrow\infty$ times.

Note: this conclusion does not rely on
1) what is the signal to be processed $d$,
2) which $F$ you use, as along as it is lowpass.

Here is some ipython code, you might play with.

import numpy as np
# generate a random low pass filter
F = np.random.randint(0,1000,(100,1) ).astype('float').ravel();
# ensures sum to 1
F = F / F.sum();
print "sum(F) = ", F.sum();
# create copy to buffer the last version of convolving F for k times
Fk = F.copy();
Fk_list = [];
for k in range( 5 ) :
    Fk = np.convolve( Fk, F, 'full' );
    Fk = Fk / Fk.sum();
    # donot worry about this line, see comments below
    Fk = np.asarray( filter( lambda x : x > 1e-6, Fk ) );
    Fk = Fk / Fk.sum();
    Fk_list.append( Fk );
# plot results
f, axarr = pyplot.subplots( 3, 2 );
# plot data
axarr[0,0].plot( F ),
axarr[0,0].set_title( 'k = 0' );
axarr[0,1].plot( Fk_list[0] ),
axarr[0,1].set_title( 'k = 1' );
axarr[1,0].plot( Fk_list[1] ),
axarr[1,0].set_title( 'k = 2' );
axarr[1,1].plot( Fk_list[2] ),
axarr[1,1].set_title( 'k = 3' );
axarr[2,0].plot( Fk_list[3] ),
axarr[2,0].set_title( 'k = 4' );
axarr[2,1].plot( Fk_list[4] ),
axarr[2,1].set_title( 'k = 5' );
pyplot.show();

the line that with filter keyword is simply to discard very small components, because I want to keep you see beautiful bell shape Gaussian. Otherwise, the results are still Gassuian, but since we might include many things outside of 4 times of standard deviation, you will see something like a peak instead.

Here are the results of the above ipython code. As you can see, the more times of convolve to $F$ self, the resulting lowpass filter is more similar to a Gaussian.

Now let me give you another visual demonstration by using a very different $F$, which only contains two nonzero values. Though this time you cannot see $F$ to be Gaussian like in 5 times of convolution, the long term behaviour namely, $k\rightarrow\infty$ does not change. It will be Gaussian like after a large number of self-convolution.

Finally, you might ask why this thing happens? Because God loves Gaussian. This conclusion is so-called the 'Central Limit Theorem' ( see details in http://en.wikipedia.org/wiki/Central_limit_theorem ). I am not gonna to explain everything for you, because I wish you could learn this on your own. But here is a hint: probability / discrete distribution, arithmetic mean, low-pass filter coefficients, linear operators, and convolution, they are all related. If you have any difficulty to understand something new, try to understand it via something you already knew with close relationships.


Does a linear phase FIR filter shifted in frequency preserve linear phase? Let's say I have a symmetric (and therefore linear phase) FIR low pass filter with real coefficients.

If I then shift this filter in some direction in frequency by multiplying its coefficients with a complex exponential (to get a band pass filter with complex coefficients), will it still be linear phase?

The generalized linear phase FIR filter has the following frequency response:

$$H(\omega) = A(\omega)~e^{j (\alpha \omega + \beta)}$$ for some constant $\alpha$ and $\beta$ and $A(\omega)$ being real. The group delay of this filter is independent of frequency: $$\tau = -\frac{d (\alpha \omega + \beta) } {d \omega} = - \alpha $$

Now if you linearly shift this filters frequency response then you get the new filter's frequency response as :

$$H_2(\omega) = H(\omega-\omega_0) = A(\omega-\omega_0)~e^{j(\alpha (\omega-\omega_0) + \beta)} = B(\omega)~e^{ j ( \alpha \omega - \alpha \omega_0 + \beta}) $$

which also shows the generalized linear phase property. The new group delay is:

$$\tau  = -\frac{d (\alpha \omega - \alpha \omega_0 + \beta) } {d \omega} = - \alpha $$

the same as the previous filter. And is also independent of frequency $\omega$.

Recall to page 295 of Discrete-Time Signal Processing 2E by A.Oppenheim for a more detailed discussion of the generalized linear phase systems.


Which is the best method for comparing the histogram obtained after BOW? Euclidian distance
Cosine similarity
Histogram Intersection
Or any other?

Consider the importance of the method in case of images. Eg, Cosine Similarity matches angle of two vectors...
What impact comparing just the angle or distance has on accuracy of BOW?
Detailed explanation is welcome!!!

Histogram Intersection (Pyramid Match Kernel) specifically is currently accepted as a better tool, but this doesn't generalize really. When you formulate your features or histogram computations differently, things are subject to change.


Inverse of Wavelet Transforms - Background and Noise Removal Main Problem: How can you inverse Wavelet Transforms
(Using the data given by signal.scipy.cwt)
I was wondering if anyone understands the scipy.signal.cwt() function well enough to use it to remove backgrounds and noise from data.
I have seen where Matlab has an inverse continuous wavelet transform function which will return the original form of the data by inputting the wavelet transform, although you can filter out the slices you don't want.
MATALAB inverse cwt funciton
Since scipy doesn't appear to have the same function, I have been trying to figure out how to get the data back in the same form, while removing the noise and background.
How do I do this?
I noticed that the wavelet has negative values, when the original data is made of Gaussian-like feautures, so I tried squaring it, but this gives me values way to large and not quite right.
I tried reading about the transforms but everything is very complicated which was difficult to parse in order to figure out what I needed to do.
Here is what I have been trying:
# Compute the wavelet transform
# I can't figure out what the width is or does?
widths = range(1,11)
# Ricker is 2nd derivative of Gaussian
# (*close* to what *most* of the features are in my data)
# (They're actually Lorentzians and Breit-Wigner-Fano lines)
cwtmatr = signal.cwt(xy['y'], signal.ricker, widths)

# Maybe we multiple by the original data? and square?
WT_to_original_data = (xy['y'] * cwtmatr)**2

And here is a fully compilable short script to show you the type of data I am trying to get and what I have etc.:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# Make some random data with peaks and noise
def make_peaks(x):
    bkg_peaks = np.array(np.zeros(len(x)))
    desired_peaks = np.array(np.zeros(len(x)))
    # Make peaks which contain the data desired
    # (Mid range/frequency peaks)
    for i in range(0,10):
        center = x[-1] * np.random.random() - x[0]
        amp = 60 * np.random.random() + 10
        width = 10 * np.random.random() + 5
        desired_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
    # Also make background peaks (not desired)
    for i in range(0,3):
        center = x[-1] * np.random.random() - x[0]
        amp = 40 * np.random.random() + 10
        width = 100 * np.random.random() + 100
        bkg_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
    return bkg_peaks, desired_peaks

x = np.array(range(0, 1000))
bkg_peaks, desired_peaks = make_peaks(x)
y_noise = np.random.normal(loc=30, scale=10, size=len(x))
y = bkg_peaks + desired_peaks + y_noise
xy = np.array( zip(x,y), dtype=[('x',float), ('y',float)])

# Compute the wavelet transform
# I can't figure out what the width is or does?
widths = range(1,11)
# Ricker is 2nd derivative of Gaussian
# (*close* to what *most* of the features are in my data)
# (They're actually Lorentzians and Breit-Wigner-Fano lines)
cwtmatr = signal.cwt(xy['y'], signal.ricker, widths)

# Maybe we multiple by the original data? and square?
WT = (xy['y'] * cwtmatr)**2

# plot the data and results
fig = plt.figure()
ax_raw_data = fig.add_subplot(4,3,1)
ax = {}
for i in range(0, 11):
    ax[i] = fig.add_subplot(4,3, i+2)

ax_desired_transformed_data = fig.add_subplot(4,3,12)

ax_raw_data.plot(xy['x'], xy['y'], 'g-')
for i in range(0,10):
    ax[i].plot(xy['x'], WT[i])

ax_desired_transformed_data.plot(xy['x'], desired_peaks, 'k-')

fig.tight_layout()
plt.show()

This script will output this image:

Where the first plot is the raw data, the middle plots are the wavelet transforms and the last plot is what I want to get out as the processed (background and noise removed) data.
Does anyone have any suggestions? Thank you so much for the help.

I ended up finding a package which provides an inverse wavelet transform function called mlpy.  The function is mlpy.wavelet.uwt.  This is the snippet which may interest people if they are trying to do noise or background removal:

# Make 2**n amount of data
new_y, bool_y = wave.pad(xy['y'])
orig_mask = np.where(bool_y==True)

# wavelet transform parameters
levels = 8
wf = 'h'
k = 2

# Remove Noise first
# Wave transform
wt = wave.uwt(new_y, wf, k, levels)
# Matrix of the difference between each wavelet level and the original data
diff_array = np.array([(wave.iuwt(wt[i:i+1], wf, k)-new_y) for i in range(len(wt))])
# Index of the level which is most similar to original data (to obtain smoothed data)
indx = np.argmin(np.sum(diff_array**2, axis=1))
# Use the wavelet levels around this region
noise_wt = wt[indx:indx+1]
# smoothed data in 2^n length
new_y = wave.iuwt(noise_wt, wf, k)

# Background Removal
error = 10000
errdiff = 100
i = -1
iter_y_dict = {0:np.copy(new_y)}
bkg_approx_dict = {0:np.array([])}
while abs(errdiff)>=1*10**-24:
    i += 1
    # Wave transform
    wt = wave.uwt(iter_y_dict[i], wf, k, levels)

    # Assume last slice is lowest frequency (background approximation)
    bkg_wt = wt[-3:-1]
    bkg_approx_dict[i] = wave.iuwt(bkg_wt, wf, k)

    # Get the error
    errdiff = error - sum(iter_y_dict[i] - bkg_approx_dict[i])**2
    error = sum(iter_y_dict[i] - bkg_approx_dict[i])**2

    # Make every peak higher than bkg_wt
    diff = (new_y - bkg_approx_dict[i])
    peak_idxs_to_remove = np.where(diff>0.)[0]
    iter_y_dict[i+1] = np.copy(new_y)
    iter_y_dict[i+1][peak_idxs_to_remove] = np.copy(bkg_approx_dict[i])[peak_idxs_to_remove]

# new data without noise and background
new_y = new_y[orig_mask]
bkg_approx = bkg_approx_dict[len(bkg_approx_dict.keys())-1][orig_mask]
new_data = diff[orig_mask]

Here a complete compilable script which can produce output:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import mlpy.wavelet as wave

# Make some random data with peaks and noise
############################################################
def gen_data():
    def make_peaks(x):
        bkg_peaks = np.array(np.zeros(len(x)))
        desired_peaks = np.array(np.zeros(len(x)))
        # Make peaks which contain the data desired
        # (Mid range/frequency peaks)
        for i in range(0,10):
            center = x[-1] * np.random.random() - x[0]
            amp = 100 * np.random.random() + 10
            width = 10 * np.random.random() + 5
            desired_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
        # Also make background peaks (not desired)
        for i in range(0,3):
            center = x[-1] * np.random.random() - x[0]
            amp = 80 * np.random.random() + 10
            width = 100 * np.random.random() + 100
            bkg_peaks += amp * np.e**(-(x-center)**2/(2*width**2))
        return bkg_peaks, desired_peaks

    # make x axis
    x = np.array(range(0, 1000))
    bkg_peaks, desired_peaks = make_peaks(x)
    avg_noise_level = 30
    std_dev_noise = 10
    size = len(x)
    scattering_noise_amp = 100
    scat_center = 100
    scat_width = 15
    scat_std_dev_noise = 100
    y_scattering_noise = np.random.normal(scattering_noise_amp, scat_std_dev_noise, size) * np.e**(-(x-scat_center)**2/(2*scat_width**2))
    y_noise = np.random.normal(avg_noise_level, std_dev_noise, size) + y_scattering_noise
    y = bkg_peaks + desired_peaks + y_noise
    xy = np.array( zip(x,y), dtype=[('x',float), ('y',float)])
    return xy
# Random data Generated
#############################################################

xy = gen_data()

# Make 2**n amount of data
new_y, bool_y = wave.pad(xy['y'])
orig_mask = np.where(bool_y==True)

# wavelet transform parameters
levels = 8
wf = 'h'
k = 2

# Remove Noise first
# Wave transform
wt = wave.uwt(new_y, wf, k, levels)
# Matrix of the difference between each wavelet level and the original data
diff_array = np.array([(wave.iuwt(wt[i:i+1], wf, k)-new_y) for i in range(len(wt))])
# Index of the level which is most similar to original data (to obtain smoothed data)
indx = np.argmin(np.sum(diff_array**2, axis=1))
# Use the wavelet levels around this region
noise_wt = wt[indx:indx+1]
# smoothed data in 2^n length
new_y = wave.iuwt(noise_wt, wf, k)

# Background Removal
error = 10000
errdiff = 100
i = -1
iter_y_dict = {0:np.copy(new_y)}
bkg_approx_dict = {0:np.array([])}
while abs(errdiff)>=1*10**-24:
    i += 1
    # Wave transform
    wt = wave.uwt(iter_y_dict[i], wf, k, levels)

    # Assume last slice is lowest frequency (background approximation)
    bkg_wt = wt[-3:-1]
    bkg_approx_dict[i] = wave.iuwt(bkg_wt, wf, k)

    # Get the error
    errdiff = error - sum(iter_y_dict[i] - bkg_approx_dict[i])**2
    error = sum(iter_y_dict[i] - bkg_approx_dict[i])**2

    # Make every peak higher than bkg_wt
    diff = (new_y - bkg_approx_dict[i])
    peak_idxs_to_remove = np.where(diff>0.)[0]
    iter_y_dict[i+1] = np.copy(new_y)
    iter_y_dict[i+1][peak_idxs_to_remove] = np.copy(bkg_approx_dict[i])[peak_idxs_to_remove]

# new data without noise and background
new_y = new_y[orig_mask]
bkg_approx = bkg_approx_dict[len(bkg_approx_dict.keys())-1][orig_mask]
new_data = diff[orig_mask]

##############################################################
# plot the data and results
fig = plt.figure()

ax_raw_data = fig.add_subplot(121)
ax_WT = fig.add_subplot(122)

ax_raw_data.plot(xy['x'], xy['y'], 'g')
for bkg in bkg_approx_dict.values():
    ax_raw_data.plot(xy['x'], bkg[orig_mask], 'k')

ax_WT.plot(xy['x'], new_data, 'y')

fig.tight_layout()
plt.show()

And here is the output I am getting now:

As you can see, there is still a problem with the background removal (it shifts to the right after each iteration), but it is a different question which I will address here.


Aquila DSP C++ Library - Wave file FFT analog frequency off by factor of 4? I'm working with the Aquila C++ DSP library. I'm computing the FFT of a wave file (16 bit depth, single channel, 44100 sample rate). I am using a window size of 16384 to calculate the FFT spectrum.

I'm having an issue converting the the digital frequency (0/N, 1/N, ... N/N) to the corresponding analog frequency. When I apply the sample rate multiplier S,

f = S*i/N; // Convert from digital frequency to analog frequency

I get analog frequencies that are 4 times greater than what I expect.

For example, I read in a pianos A1 note. I expect this to be ~55 Hz. Instead I'm getting ~220 Hz. Refer to the photo below.

Now I'm first wondering if I'm scaling my the wrong number. However, I'm pretty sure it's correct. So now I'm curious if somehow the time domain data I'm reading is compressed by a factor of 4, such that the frequencies appear to be 4 times greater. Does anyone know if there is a common .wav file 'gotch-ya' when reading from them?

I've posted the relevant part of my code if anyone is interested.

/* sample window */
int window_size = 16384;

/* Calculate the FFT */
std::shared_ptr<Aquila::Fft> p_fft_interface = Aquila::FftFactory::getFft(window_size);  // This returns a shared pointer to an FFT calculation object.
auto spectrum = p_fft_interface->fft(wave_object.toArray());

QVector<double> x(window_size);
QVector<double> y(window_size);

/* Prepare to plot & convert to analog frequency */
double max_value = 0;
for(int i = 0; i < window_size; i++)
{
    x[i] = i*(sample_freq/window_size);
    y[i] = abs(spectrum[i]);

    if(abs(spectrum[i]) > max_value)
    {
        max_value = abs(spectrum[i]);
    }
}

Not familiar with Aquila in particular, and I don't see how you're mapping the output spectrum samples to an X-axis that ends at 800 Hz,  but remember that a real-only FFT mirrors its output data in the second half of the buffer.  You may need to run your loop from 0 to window_size/2 if Aquila's FFT is implemented the usual way.

That could account for a factor of 2x.


Can I compute convolution through the overlap and add method without the IFFT? So, I'm (I think) aware of how the overlap and add algorithm for linear convolution works, but my question is that, suppose I have a FFT-ed set of sequences that belong to a large sequence.

Can I add the overlapping parts before taking the IFFT of each sequence?
Can I add the overlapping parts and then perform another convolution without taking another set of IFFTs and FFTs?
Let's say that I have a function which multiplies by 3 if the number is 0 in the time domain, but multiplies by 6 if it isn't. I can perform this with only one IFFT as opposed to an IFFT+ a FFT by taking the IFFT and then multiplying in the frequency domain by the required number. Can I do a convolution through overlap and add, apply the aforementioned function through the "one IFFT only technique" and then do another convolution through the overlap and add method without having to do another set of FFTs for the next convolution through the overlap and add?

your FIR, $h[n]$ can be converted, once and for all, into $H[k]$ the spectrum with which you will multiply in the frequency domain.  but for every frame of $x[n]$, you must zero pad and FFT that into $X[k]$, multiply that $X[k]$ by $H[k]$ (which now becomes $Y[k]$) and you must iFFT that back into $y[n]$ and overlap add with the previous frame of $y[n]$.


Has this 2D filter for enhancing circular dots in images a name? I came across this 2D filter for enhancing circular dots in images, for example to enhance a dot with a diameter of 5 pixel, the filter is:

$\frac{1}{336}\left( \begin{array}{ccccccc}
0 & 0 & −21 & −21 & −21 & 0 & 0\\
0 & −21 & 16 & 16 & 16 & −21 & 0\\
−21 & 16 & 16 & 16 & 16 & 16 & −21\\
−21 & 16 & 16 & 16 & 16 & 16 & −21\\
−21 & 16 & 16 & 16 & 16 & 16 & −21\\
0 & −21 & 16 & 16 & 16 & −21 & 0\\
0 & 0 & −21 & −21 & −21 & 0 & 0 \end{array} \right)$

Has this filter a name? Like Sobel or Scharr.

What is the rationale behind its coefficients?

This is a matched filter, where the shape of the filter coincides with the shape of the signal to be detected.

In this case, whenever the filter convolution overlaps a circular dot in the image, a maximum will be present in the output image. Thus, searching for the single maximum pixels in the output image will give the center position of dots in the original image.

Note that your linked source mentions that it is a matched filter.


Detecting text lines in image by means of Hough Transform In my Android App,I have to detect text lines in images.I want to use Hough Transform to detect lines.You can see an image below:

I got edges in image by canny filter,then I applied Hough Transform on canny image:

Size size = canny.size();
Mat lines = new Mat();
int lineGap = 20;
double sum1 = 0;
Imgproc.HoughLinesP(canny, lines, 1, Math.PI / 180, 100, size.width / 2.f, lineGap);

Result for above image, looks like this:

You can see three text lines are detected.So I added line's thickness by applying dilate on canny image:

Mat dilate = new Mat();
Mat kernel = Imgproc.getStructuringElement(Imgproc.MORPH_RECT, new Size(10,10));
Imgproc.dilate(canny, dilate, kernel);
int lineGap = 20;
double sum1 = 0;
Imgproc.HoughLinesP(dilate, lines, 1, Math.PI / 180, 100, size.width / 2.f, lineGap);

I drew detected lines on empty black background.Result is shown below:

You can see number of detected text lines is decreased,although I had added thickness of lines and did not change amount of minLineLength or maxLineGap,also when I apply dilation no white pixel remove from image!

Again I added thickness of lines by doing dilation with kernel of size(30,30):

Mat dilate = new Mat();
Mat kernel = Imgproc.getStructuringElement(Imgproc.MORPH_RECT, new Size(30,30));
Imgproc.dilate(canny, dilate, kernel);
int lineGap = 20;
double sum1 = 0;
Imgproc.HoughLinesP(dilate, lines, 1, Math.PI / 180, 100, size.width / 2.f, lineGap);

Result is:

Now my question is:
Why in each level when I add number of white pixels (by dilation) and do not change minLineLength or maxLineGap,Hough Transform do not detect some lines that were detected in previous level with less white pixels?

OpenCV now has a text detection module included. You might want to take a look at it:

Detector: http://docs.opencv.org/master/modules/objdetect/doc/erfilter.html
Recognizer: http://docs.opencv.org/trunk/modules/text/doc/ocr.html


Why root raised cosine filter can eliminate intersymbol interference (ISI) ? My question is kinda silly, I know it. But I cannot understand why root raised cosine filter or raised cosine filter can eliminate ISI.
I googled before, but those explanation is not easy for me to understand. Maybe someone can help me this with some figures? Thank you so much for you help in advance.

In addition to the frequency response that AlexTP showed, I think the time domain view of the impulse response gives good insight into the Zero-ISI operation of a Nyquist filter. Below is shown the result of transmitting three symbols, than each began as a single impulse. Each plot below is the "impulse response" for each symbol. Note the blue symbol was transmitted first and is maximum at time index 0 on the graph (but began much earlier than this in time). This position where the impulse response is maximum is our ideal sampling location to determine what was transmitted. One sample after the original impulse for the blue symbol, came an impulse for the red symbol. Notice at the ideal sampling location for red (where it is maximum), the response from all the other symbols is zero (meaning zero ISI).

The pulses shown above have zero ISI at the sampling locations, but extend to infinity to achieve the maximimum conservation of bandwidth (1/(2T)). Of course, this is not feasible, so the raised cosine pulse shape applies a window to truncate the perfect Sinc function, resulting in a trade of bandwidth versus filter complexity, as well as trades with other characteristics related to performance for timing recovery.

The filter is often split between transmitter and receiver as a "root-raised cosine filter" (RRC) on each side. Note that a RRC filter does NOT have zero-ISI until it is paired with the second RRC filter to form in cascade a raised cosine filter.
Below is an example eye diagram, showing the trajectory over 1 symbol for multiple pulses (each trajectory is based on the history of all previous symbols within the depth of the RRC filter and results in a different path when going from 0 to 1 or 0 to 0 etc).  At time 0 is the ideal sampling instant to decide between 0 and 1, notice with the RRC filter that there is significant ISI, which is eliminated after passing through the second RRC filter in the receiver (assuming no other multipath effects have occurred).  Also notice for the filter with alpha = 0.25 that there is a significant spread in the zero crossing locations, which would add more noise in timing recovery operations that are sensitive to zero crossings. Compare this to the eye diagram with alpha= 1 and notice the significant decrease in zero-crossing jitter, as well as a wider "eye-opening" which also means less sensitivity to jitter or other errors in the sampling clock location, but this benefit comes at the expense of less spectral efficiency.

What is very interesting is this alternate pulse shaping Nyquist filter in this paper http://eon.sdsu.edu/~seshagir/SDR05.pdf that fred harris now refers to as the "harris-Moerder Nyquist Pulse". Using this as the pulse shaping filter results in an order of magnitude lower EVM for the same number of taps. This isn't commonly used even today (15 years later) to my knowledge since the root-raised cosine filter is baked into our specifications but certainly is very compelling to be aware of for consideration in future radios.


Manipulating RGB color by copying channel value I have a picture and trying to swap the RGB channel into RRB, it's like I need to change all G pixels value into R. I did that by copying the R pixel value to G and store it into a new image file (Mat).
Here's my code :

ori.at<Vec3b>(i, j)[0] = image.at<Vec3b>(i, j)[0]; // B-->B
ori.at<Vec3b>(i, j)[1] = image.at<Vec3b>(i, j)[2]; // G-->R
ori.at<Vec3b>(i, j)[2] = image.at<Vec3b>(i, j)[2]; // R-->R

however the result is a little bit confusing me,

Could someone give me a brief explanation how is that the Red and Yellow pixels in the wall end up have the almost same result color (aquamarine) ? I checked the pixel value at some points and the value is almost the same.

What I expected was, if I copied R pixels --> G, part with original R values shouldn't be affected, right ? And because Yellow is a mixture of Red + Green, if I try to copy Red pixels into Green, shouldn't it be Red + Red then ?

[UPDATE]

I have checked sample color values (RGB) for Red wall and Yellow wall.

Red wall :              Yellow wall :
R --> 200 195 185       R --> 237 226 220
G -->   2   1   1       G --> 197 179 168
B -->   1   1   0       B -->   4  29  23

Copying the value R --> G :
Red wall :              Yellow wall :
R --> 200 195 185       R --> 237 226 220
G --> 200 195 185       G --> 237 226 220
B -->   1   1   0       B -->   4  29  23

will result in the Red wall and Yellow wall have almost the same color. But when I checked the result color, it's not aquamarine. Did I make wrong calculation ?

I missed the conversion part.
I converted the result image from BGR into RGB using this command :

cvtColor(image, image, COLOR_BGR2RGB);

Then I exported the result image into image file.

imwrite("image.png", image);

Turned out the export command will automatically convert the BGR color-space into the RGB color-space. This is why the color result is off because I did 2x conversion.
So, I don't need the conversion command.


A question about sampling theory Let`s say I have a signal

$m(t)=\cos(4\pi t) + \cos(6\pi t)$

so we can say the signal is containing frequencies $f_1= 2\mathrm{Hz}$ and $f_2= 3\mathrm{Hz}$. The cut-off frequency of the low pass filter is equal to $f_c= 3.5\mathrm{Hz}$, and the sampling frequency is $5\mathrm{Hz}$.

Is there an overlap (aliasing error) when receiving the original signal (After LPF)?

Given two signals $x_1(t) = \cos(2\pi 2t)$ and $x_2(t) = \cos(2\pi 3t)$ which have frequencies of 2 and 3 Hz respectiveley, if you sample them uniformly with a sampling rate of 5 samples per second (or a sampling period $T$ of $\frac {1}{5}$ = 0.2 seconds) you will get the following two discrete time signals $$x_1[n] = x_1(nT) =  \cos(4\pi \frac {n}{5} )$$ and $$x_2[n] = x_2(nT) = \cos(6\pi \frac {n}{5})$$

It can be easily seen with simple trigonometry, ( by using $\cos(x) = \cos(2\pi-x)$) ,  that $x_1[n] = x_2[n]$ which is a manifestation of aliasing that occured on the 3hz signal, therefore, the 3Hz signal seems like 2Hz.


How to find a 'square root filter' such that $d*\overline{d}=f$ given $f$? I have a desired frequency response $F(\omega)$. I know how to create a FIR filter (say, with linear phase) $f(n)$ with this frequency response.

But how can I get a filter $d(n)$ such that $D(\omega)D^*(\omega)=|D(\omega)|^2=F(\omega)$, or equivalently in the time domain, $(d*\overline d)(n)=f(n)$, i.e. the 'square root' of $f$? (here $\overline d$ is the filter $d$ with reversed order of coefficients)

I would like to get a real filter with compact support, preferably using MATLAB.

Thanks.

If you can design an FIR filter with arbitrary magnitude response $F(\omega)$, then you can easily design the filters that you're looking for. Just break the design into two pieces: magnitude and phase.

If the desired magnitude response of the cascade is $|F(\omega)| = |D(\omega)||D^*(\omega)| = |D(\omega)|^2$, then, as was pointed out already, you can simply take the square root of the desired overall magnitude response to yield the magnitude response of your square-root filter $D(\omega$).

With the magnitude response defined, then it's up to you to pick what kind of phase response you want to give $D(\omega)$. The phase responses of $D(\omega)$ and $D^*(\omega)$ will be conjugates of one another, so they will cancel when the two square-root filters are cascaded. You would typically pick the phase response to aid in implementation or to give you some other nice property along the way.

One common choice, as you pointed out, is using linear-phase FIR filters. So, you can pick a filter length, set the filter group delay to $\frac{N}{2}$ (where $N$ is the filter order), and use the resulting phase delay curve as the desired phase response in your square-root filter design. You should end up with a filter with a symmetric impulse response (which can provide implementation advantages due to the redundancy in the structure of the taps).

One other thing I should note is that the question as posed will yield a non-causal structure in one of the subfilters. Take the case where $D(\omega)$ has linear phase (so a constant group delay). Its conjugate will have the opposite phase response, so if $D(\omega)$ corresponds to a delay of $N$ samples, then $D^*(\omega)$ will have a phase response that corresponds to an advance of $N$ samples, so it will be non-causal.

What the OP was probably looking for was a structure where $D^2(\omega) = F(\omega)$ instead. This would be similar to the root-raised cosine filter structure used in digital communications, for example.


Pa^2/Hz to dB/Hz conversion I have PSD data from an computational aeroacoustics analysis, which is in the units of Pa^2/Hz. I want to convert this to dB/hz; would it be correct if I use PSD[in dB/Hz] = 10log10(PSD[in Pa^2/Hz]) or is there another way of doing this that I am not aware of?

This depends a bit on what you are after. As Laurent Duval has pointed out, you need can only take the logarithm of a unitless quantity so you need to do some form of normalization

If you want Sound Pressure Level (dbSPL) you would normalize to the standard reference pressure of 20 uPa (20e-6 Pa)
You could also divide by the free field acoustic impedance of air, which would give you the intensity in W/m^2.
The most common units for a PSD would be W/Hz. To calculate that you'd have to integrate over the entire volume flow and so you need to make an assumption for the shape of your sound wave. That's pretty easy for spherical and plane waves, but difficult for more complex shapes.


What part of complex number of inverse discrete Fourier transform? Ok, so we have an image that is a Fourier inverse of the original picture. We want to get the original picture back. We use Matlab to get that job done. We import the image and then we invert it with the help of ifft(), this gives us a matrix with complex numbers. But to get the original picture we need to do some operation on the complex numbers to get it. But what is that operations. I tried the magnitude, real and imaginary part but this doesn't create the picture we want.

To apply IFFT you need back the signal do complex numbers, you need use magnitude and phase information to rebuild correctly.

The real part is = magnitude * cos(phase)

The imaginary part is = magnitude * sin(phase)

You can use square roots of −1 (sqrt(-1)) to get Imaginary unit.

Now multiply imaginary unit with imaginary part and sum with real part, OK now are you done to apply IFFT !

At the end I apply a mat2gray function to convert the matrix to the intensity!

here how it is really done in matlab:

x=imread('C:\Users\Eder\Pictures\download.jpg');
figure(1);imshow(x);
%Make FFT
y=fft(x);
%Amplitude of the FFT
mx=abs(y);
%get Phase Information
ma=angle(y);
%back the signal to complex
y2= mx .* ( cos(ma) + sqrt(-1) * sin(ma) );
%Apply Inverse FFT
x2=real(ifft(y2));
result=mat2gray(x2);
figure(2);imshow(result);


How is "Joint class histogram" algorithm deploy using Matlab? Image registration and image retrieval depend on image similarity measure to find mutual information between 2 images. Using histogram, joint histogram, and joint class histogram we can find the degree of similarity between those images. Hopefully, I understand and deploy both first techniques.Here, I found joint histogram algorithm and deploy using Matlab.

For example:
this paper was mention this histogram and used as image similarity measurement, but I couldn't found how the authors calculate it.

However, I search for many days to find good paper, book or article regarding joint class histogram but their authors only mention this technique without enough detail. I think this technique still infant and doesn't have enough study yet.

I think joint class histogram algorithm will the best answer for my question.

I get the answer from another source and I deploy like this

% load floating and fix image

%% image segmentation using k-mean clustring.
idx_fix = kmeans(double(fix(:)),4);
idx_float = kmeans(double(float(:)),4);

%% display image segmentation
mat1 = reshape(idx_fix,[],500);%besed on your image size
mat2  = reshape(idx_float,[],500);
imshow(mat1/max(mat1(:))); title('Fix image class');
figure;
imshow(mat2/max(mat2(:))); title('Float image class');

%% Joint class histogram
joint = zeros(5,5);%besed on number of class
for i=1:length(idx_fix)
   c1 = idx_fix(i) + 1;%offset matrix by 1 because index start by 1
   c2 = idx_float(i) + 1;
   joint(c1, c2) =  joint(c1,c2) + 1;
end

If the image aligned correctly the joint class histogram result as a diagonal line.

Thanks, Walter Roberson


How to Compute Image Gradients? How to compute Image gradients?

In X Direction - Kernel : -1 0 1

In Y Direction - Kernel :
-1
0
1

Let's assume my input image is: unsigned char* imgImage;

Could any one please share code. I am newbie to image processing.
Is it possible to do in OpenCV with this same kernel?

Thanks a lot in advance!

You can create your own custom kernel filter using something similar to this example.

if you just want to find image gradients there are other options such as sobel and laplace

If your aim is edge detection, I find canny is best for this in most cases.


Why we need distances in Ideal LPF? Here I have a MATLAB code for Ideal Low pass filtering an Image. In the below code they have used D=sqrt(U.^2+V.^2); which finds some distances. Then this distance matrix D is compared with the cut-off frequency to create a filter function H ( H=double(D<=P); ) that has be multiplied with the fourier transformed image F. i.e, G=H.*F; to get the output image.

My questions: What is this distance ? Why we need this distance matrix D ? What information this gives to us? How it helps in Low Pass Filter an Image? (I understood that it helps to create the filter function H) but Why we need the distances to create the filter function? I need to it understand clearly. This may be a simple question but could anyone please explain me this with little patience.

Thank you for reading my question.

Reference: MATHWORKS:IDEAL LOW PASS FILTER

%IDEAL LOW-PASS FILTER

    function idealfilter(X,P) % X is the input image and P is the cut-off freq
    f=imread(X);  % reading an image X
    [M,N]=size(f); % Saving the the rows of X in M and columns in N
    F=fft2(double(f)); % Taking Fourier transform to the input image

    u=0:(M-1);
    v=0:(N-1);
    idx=find(u>M/2);
    u(idx)=u(idx)-M;
    idy=find(v>N/2);
    v(idy)=v(idy)-N;
    [V,U]=meshgrid(v,u);
    D=sqrt(U.^2+V.^2);   % finding out the distance

   H=double(D<=P);       % Comparing with the cut-off frequency
    G=H.*F;               % Multiplying the Fourier transformed image with H
    g=real(ifft2(double(G))); % Inverse Fourier transform
    imshow(f),figure,imshow(g,[ ]); % Displaying input and output image
    end

In the Fourier domain, the frequency increases with the distance from origo. So an ideal low-pass filter in the Fourier domain is shaped like a disk centered at the origo.

Here they work with only one quadrant of the Fourier domain. Therefore, they create a pie-shaped mask H for the quadrant by computing the Euclidean distance from the quadrant corner and store the distance values in D. The mask H is then chosen as all elements in D with values less than the parameter P (which is the radius of the disk).

So the answer to your question(s) is that the distance computation helps the method compute the filter shape. The key is that a constant Euclidean distance from a point defines the shape of a circle, and all elements inside that circle form the disk that is used as filter mask.

So, naturally, after filtering, only the low-frequency components of the image are left in the Fourier domain. Then the low-pass filtered image is retrieved by inverse transforming the filtered image back from the Fourier domain to the spatial domain.

EDIT

Here is an illustrated example to clarify the answer.

A typical way of illustrating filtering in the Fourier domain is to show the input image.

Then we show the Fourier transform of the image.

Here the origo of the Fourier space is in the center. Since the frequencies increase with the distance from origo, the ideal filter is a disk with its center at the origo. Here we see the disk.

The disk radius is the input parameter P. We use the disk to filter the image by multiplying the filter with the transformed image. We then get the filtered transform.

To get the filtered image we take the inverse transform of the filtered result.

So far so good. As I said, this is usually how low-pass filtering in the Fourier domain is explained.

Now what I was trying to explain when talking about quadrants of the Fourier domain, which got a bit confusing, is that the result from the Fourier transform in the MATLAB code provided does not have origo centered as in the second image in the example above (the image showing the unfiltered Fourier transform).

Instead the result from the Fourier transform F=fft2(double(f)) looks like this (we refer to it as the transform image).

This means that since the filter is a disk with its center in origo, the disk must be placed in the corner of the transform image. Since the Fourier transform is cyclic, the filter should wrap around to cover all four corners of the transform image.

Here is where the distance matrix D comes in. It is used to compute the filter mask above. D is computed by taking the Euclidean distance from each corner (using a neat trick with meshgrid). So Dis the Euclidean distance map looking like this,

where blue values represent low distances and red values represent high.

% Top left corner of D
0.00000   1.00000   2.00000   3.00000   4.00000   5.00000
1.00000   1.41421   2.23607   3.16228   4.12311   5.09902
2.00000   2.23607   2.82843   3.60555   4.47214   5.38516
3.00000   3.16228   3.60555   4.24264   5.00000   5.83095
4.00000   4.12311   4.47214   5.00000   5.65685   6.40312
5.00000   5.09902   5.38516   5.83095   6.40312   7.07107

The filter mask can therefore simply be acquired by selecting the elements in D that have a distance value less than or equal to P.

The rest of the code works as in the example described above. That is, the filter mask is multiplied by the transform image, and the result is inverse transformed to get the low-pass filtered image.

I hope this clarifies what role the distance computation has for creating the filter.


On the stability and causality of a discrete system On  MIT's open course a simple exercise with two questions is given. On the first part, they question about the properties of the following discrete system:

$$
y[n]=x[n]+0.5y[n−1]−2y[n−2]
$$

The answer states that none of the following properties can be determined, as no initial conditions were given: stability, causality, linearity, time (in)variance.

I would be inclined to say that the system is at LTI, as the coefficients are constant and the output is simply defined by a recursive relationship to itself and the input. Both poles are inside the unit circle and the system should be stable if it is causal, as the ROC includes the unit circle.

My first question is therefore: could you come with a(n) example(s) of initial conditions (or any argument for that matter) so that the system is: variant in time OR unstable OR not causal OR non-linear? The only conditions that I can think of that could make it behave strangely would be extreme ones such as +/- infinity. I assumed this was not the intent of the initial question.

Next, the following transfer function is given (unrelated to question above):

$$
H(e^{j\omega }) = \frac{1 - 2e^{-j\omega }}{(1 - \frac{3}{4}e^{-j\omega })(1 - 3e^{-j\omega })}
$$

Here it is also asked about its properties. On the final answer they state that, given that a transfer function is given, one can imply that the system is linear and time invariant. So far so good. Next, however, it states that the system is both stable and not causal, with the following argument:

...Further, from the existence of a frequency response you may conclude that the frequency response converges on the unit circle, i.e., that the system is stable. Given that the system is stable the region of convergence is also implied resulting in a two-sided sequence, meaning that the impulse response is not causal.

Hence my second question:

How was this concluded from this transfer function, given that no ROC constraint is there?

Thanks in advance.

To your first question, a linear system is causal if we have initial conditions of rest that is zero. This is a sufficient and necceassy condition. This property doesn't  require time invariance

Stability is a consequence of absolute summability of system response. It doesn't depend on the initial condition (of course have to be finite to make any sense).

To your second question the transfer function has a pole at z=3 and it is given that the system is stable. That means the ROC has to include the unit circle. That means that the ROC is to the left of the pole at z = 3, since ROC is towards the left of outermost pole hence system is not causal


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

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

  where $n$ and $m$ are integers.

Using the definition I get

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

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

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

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

or, alternatively

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

and you will need to use one of those.

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

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

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

for autocorrelation of a purely real signal

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


How to calculate the steady state response $y_{ss}(t)$ of a LTI system given the Laplace transform $Y(s)$? I am given the Laplace transform of the output of a LTI system: $$Y(s) = \frac{1}{s((s+2)^2+1)}$$ Asked is what the steady state response $y_{ss}(t)$ would be. I think that $y_{ss}(t) = \lim_{t\to\infty} y(t)$, since after waiting infinit long, the system should be in steady state. (Right?)

I thought to use the final value theorem:

$$\lim_{t\to\infty}y(t)=\lim_{s\to 0}sY(s)$$
This gives $$\lim_{s\to 0} sY(s)=\lim_{s\to 0}\frac{1}{(s+2)^2+1} = \frac{1}{5}.$$

This is different from $\frac{1}{10}$. When I let a computer algebra system calculate $\mathscr{L}^{-1}[Y(s)]
\bigg{|}_{t=\infty}$ I get $\frac{1}{10}$.
(I'm using wxMaxima and used limit(ilt(1/(s*(s^2 + 2*s + 10)), s, t), t, inf);.)

What am I doing wrong? Thanks in advance.

For these calculations, it is better to give the Wolfram Alpha answers:

inverse Laplace transform 1/(s*(s^2 + 4*s + 5))

Which gives the correct expression, consistent with the Final Value expression of 1/5:

$$\frac{1}{5} - \frac{1}{5} e^{-2 t} (2 sin(t) + cos(t))$$

This is very different from the complex (phasor) expression, which could inexplicably have a 1/10, but when evaluated it is still 1/5:

$$\frac{1}{10} i (e^{(-2 - i) t} ((2 + i) e^{2 i t} + (-2 + i)) - 2 i)$$


In signal averaging, will the noise be reduced if the noise is not Gaussian? My understanding is that signal averaging works this way. If we take the average of the waveform, we get the average of the signal component (which is the magnitude of the signal) and the average of the noise component (which becomes closer to zero as number of samples are increased).

However, this is based on the assumption that the noise has zero mean and constant variance (i.e. white or Gaussian noise). What if the noise is not white? For example, what if the noise decreases as frequency increases? Will the noise still be averaged out by signal averaging?

Consider two zero mean signals $x(n)$ and $y(n)$. Let $z(n)=x(n)+y(n)$. In the following $E()$ denotes the expectation.
$$E(z) = E(x+y)= E(x)+E(y)= 0$$
I've dropped the time index $n$ for convenience. So $z(n)$ has zero mean. Now let's look at the variance. Since $z(n)$ is zero mean - I'll just use the correlation.
$$E[z^2] =E[(x+y)(x+y)]=E[x^2+2xy+y^2]=E(x^2)+2E(xy)+E(y^2) $$

If the correlation between $x(n)$ and $y(n)$ is zero, then $E(xy)=0$, then the variances add. Usually this is done by stating the $x(n)$ and $(y)$ are independent which implies uncorrelated.

In the general case $E(xy)\neq 0$, so the resulting variance could increase or decrease depending on the sign of the correlation.

Note - stating that $x(n)$ and $y(n)$ are white can be a little unclear. Usually saying $x(n)$ is white implies it is uncorrelated with itself - it does not state how it is correlated to another noise process.

The above analysis is not dependent on the type of noise (Gaussian or not) and it does not make any assumptions on the autocorrelation functions of $x(n)$ and $y(n)$.


How to calculate DTFT of cosine function divided by n I'm having a hard time to calculate the next function, and I don't really know Matlab good enough to calculate it there.
Help would be appreciated:
$$h[n]=\frac{A_1 \cos[\theta_1(n-N/2)]}{n-N/2}$$

The given impulse response is a scaled version of the impulse response of a highpass Hilbert transformer with delay $\tau=N/2$ ($N$ odd) with frequency response
$$H(e^{j\omega})=\begin{cases}-j\,\textrm{sgn}(\omega)\,e^{-j\omega\tau},&\omega_c<\omega<\pi\\j\,\textrm{sgn}(\omega)\,e^{-j\omega\tau},&-\pi<\omega<-\omega_c\\0,&-\omega_c<\omega<\omega_c\end{cases}\tag{1}$$
This can be shown as follows:
$$\begin{align}h[n]&=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\omega})e^{jn\omega }d\omega\\&=\frac{1}{2\pi}\left[-j\int_{\omega_c}^{\pi}e^{j\omega(n-\tau)}d\omega+j\int_{-\pi}^{-\omega_c}e^{j\omega(n-\tau)}d\omega\right]\\&=\frac{1}{2\pi (n-\tau)}\left[-e^{j\pi(n-\tau)}+e^{j\omega_c(n-\tau)}+e^{-j\omega_c(n-\tau)}-e^{-j\pi(n-\tau)}\right]\\&=\frac{\cos[\omega_c(n-\tau)]-\cos[\pi(n-\tau)]}{\pi(n-\tau)}\tag{2}\end{align}$$
Since for $\tau=(2k+1)/2$, $k\in\mathbb{Z}$, we have $\cos[\pi(n-\tau)]=0$ we arrive at the following inverse DTFT of $(1)$:
$$h[n]=\frac{\cos[\omega_c(n-\tau)]}{\pi(n-\tau)}\tag{3}$$


Interpreting N in DFT as the Number of Points vs. Number of Intervals The "N" is DFT is understood to be the number of data points in a given sequence or in other words the length of the sequence. We recently have had discussions here Indexing in DFT (from an old paper) and someone's old question How do I measure the time duration of a finite-length discrete sequence?. One of the popular symbolic versions of DFT is
$$
X(j)=\sum_{k=0}^{N-1} x(k) \exp \left(-i 2 \pi\left(\frac{j}{N}\right) k\right)
$$
Suppose someone gives us the sequence consisting of N=11 points without telling us the total time or the sampling rate. If we apply DFT on it in MATLAB, the output is 11 points

Notational problems start when we wish to determine the time interval $\Delta$t and the frequency step $\frac{1}{N\Delta t}$ when the $\Delta$t is revealed.
a) If we wish to determine the sampling frequency it is (N-1) points collected in 1 seconds. The last point belongs to the next cycle. The correct sampling rate is 10 Hz not 11.
b) The frequency resolution is shown to be $\frac{1}{N\Delta t}$. In order to get the correct  frequency step, we to have to put $\frac{1}{10(0.1)}$  NOT $\frac{1}{11(0.1)}$.
Therefore, the frequency resolution, if we accept that N=11,
$$\frac{1}{(N-1)\Delta t}$$ but nobody shows this formula.
It seems that we are using N in two ways

N as the length of the sequence
N when we have to determine the frequency step where actually it is one less than N to get the correct result.

A respected colleague says that N should be interpreted as the no of intervals not the number of points. This is inconsistent with defining N as the number of points in the sequence. Unfortunately, I cannot find any reference which says that N is the number of intervals.
How can we make this consistent?
Thanks.

The key is to understand what the DFT says, vs what we seek. Consider a cosine, where we change $f \text[Hz]$, $N$, and $t$ and observe the effect on DFT:

[1]: DFT "sees" 1 cycle in the "analysis frame" (i.e. what we feed it), so nonzero at $k=1$, as expected.
[1] to [2]: we double our time duration without changing $f$ or $N$; DFT sees this as two cycles spanning the analysis frame, so $k=2$.
[1] to [3]: we double the physical frequency without changing $N$ or $t$; DFT sees this as two cycles spanning the frame.
[1] to [4]: we double the number of samples, $N$, without changing $f$ or $t$; DFT sees this as still 1 cycle spanning the frame; bin location, $k$, remains unchanged, but (unnormalized) correlation strength doubles (not relevant here).
[4] to [5]: now we also double duration, yielding two cycles in analysis frame.

You should be seeing a pattern. Without reading further, try to establish a relationship between $k$, $N$, $t$, and $f$. Hint: units.

Here's the deal: DFT has no idea what Hz, or physical frequency, is. All it knows is samples and cycles spanning the analysis frame. In [1], the "DFT frequency" is
$$
f_{\text{DFT}} = \frac{k}{N} = \frac{1 \text{ cycle}}{10 \text{ samples}} = .1 \left[ \frac{\text{cycles}}{\text{samples}} \right]
$$
[2] = [3] = 2 cycles / 10 samples, [4] = 1 cycle / 20 samples, [5]  = 2 cycles / 20 samples. Now let's take what we know about physical frequency, $f_p$ and DFT frequency, and relate them. In [2], DFT says $k=2$, but we know $f_p = 1$. It could also be (not shown in any [1]-[5]) $f_p=2$ and $k=1$. How to convert?
The unifying relation is:

 \begin{align} f_p \left[ \frac{\text{cycles}}{\text{second}} \right] & = \left( f_{\text{DFT}} \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right) \end{align}

So, for [2]:
$$
f_{\text{DFT}} \cdot f_s =  \left( \frac{2 \text{ cycles}}{10 \text{ samples}} \right) \cdot \left( \frac{10 \text{ samples}}{2 \text{ seconds}} \right) = 1 \left[ \frac{\text{cycles}}{\text{second}} \right] = 1\ \text{Hz} = f_p
$$

But how to determine $f_s$?
By a simple definition, it's the reciprocal of sampling period, $\Delta t$, making everything above consistent. However, one must ask, if "sampling frequency" is defined as "# of samples / total duration", and "total duration" of
$$
[0, .1, .2, .3, .4, .5, .6, .7, .8, .9]\ \text{sec}
$$
is clearly $0.9\ \text{sec}$, then isn't $f_s$ actually $.9 / 10 = 0.9\ \text{Hz}$? No; the duration is actually 1 sec. Here's why: $0.9\ \text{sec}$ here is actually the duration of something else entirely. Namely, "what is the duration of the signal?" can ask two things:

For how long have we been sampling?
How much time's worth of information is contained in our signal?

The answer to former is $0.9\ \text{sec}$, but to latter is $1\ \text{sec}$. Former computed via $(N-1)\Delta t$, latter via $N \Delta $, and if we insist on 0.9 for #2, we're saying that one sample holds no time-representative information, which implies all signals have duration of zero.
I clarify this with examples here. In a nutshell, the end-goal is description of information, not of process used to obtain it.

So what is the frequency resolution (rather DFT bin spacing)${}^{1}$?
It's defined as spacing between DFT bins, $df$; the answer depends on units of choice. For Hertz, per all of the above,
$$
df_p = \frac{1}{N \Delta t} \tag{1}
$$
meaning, for [2], $k=1$ corresponds to $f_p = 0.5\ \text{Hz}$, $k=2$ corresponds to $f_p = 1\ \text{Hz}$, and so on. Alternatively, if you insist on defining duration via $(N-1)$, then it'll be per $(N-1)$ in Hertz, but not in DFT frequencies; latter is unambiguous:
$$
k=1 \rightarrow \frac{1\ \text{cycle}}{N\ \text{samples}} = \frac{1}{N} \left[ \frac{\text{cycles}}{\text{sample}} \right]
$$
You can again convert between DFT frequency resolution and physical; taking [2], $\text{Duration} / N = 2 \text{ sec} / 10 = .2 \text{ sec}$, so bin spacing is
$$
df_p = \frac{1}{N \Delta t} = .5\ \text{Hz}
$$
The spacing can change, but it's by redefining $\Delta t$ rather than changing $N$ to $(N - 1)$ in $(1)$. Suppose we say $\text{Duration} = 1.8\ \text{sec}$; then, $\Delta t = 0.18\ \text{sec}$, and
$$
df_p = \frac{1}{N \Delta t} = 0.\bar{5}\ \text{Hz}
$$
So in [2], $k=1$ corresponds to $0.56\ \text{Hz}$, and $k=2$ corresponds to $1.1\ \text{Hz}$, which agrees with completing 2 cycles in 1.8 secs = 1.1 Hz.

Your specific example: $N=11$, $\Delta t = 0.1\ \text{sec}$:
$$
df_p = \frac{1}{N \Delta t} = \frac{1}{11 \cdot 0.1\ \text{sec}} = 0.909\ \text{Hz}
$$
So $k=1$ corresponds to $0.909\ \text{Hz}$, not to 1Hz, since you've included a sample from the next cycle in the analysis frame.

1: NOTE: $df$ is DFT bin spacing, not "frequency resolution". DFT has perfect frequency resolution and no time resolution. But if you define it as discrimination of continuous-time frequencies, then the resolution and bin spacing are inversely-related (lesser spacing -> more bins -> more granular resolution). This is its own topic so I'll avoid clarifying in detail, feel free to open new q.


Is Entropy Coding in Wavelet equivalent to Huffman Coding in DCT How exactly we should describe the process of Entropy Coding in wavelets ? can't it be simply called Coding ?

Here is the block Diagram of Image compression using wavelet,

I am confused about the "Entropy Coding" part. Can i simply use Huffman coding to code bits ?

The goal of entropy coding is to code your data such that the length of the code associated with a symbol is close to $-\log_2 p$, where $p$ is the probability of the symbol (Shannon's limit). This can be achieved through different techniques: Huffman coding is one of them, but arithmetic coding is another option - or if the distribution of the coefficients can be known in advance, a pre-computed code can be specified. The choice of an entropy coding method is motivated by performance/complexity trade-offs, and by the known (or unknown) distribution of the coefficients. If two transforms generate the same distribution of coefficients, the same entropy coding methods can be used for both.

Since the diagram from your question hints that your application is JPEG2000 compression, I suggest you to look into the JPEG2000 specifications to see which entropy coding method they have actually specified.


3CCD or monochrome with Bayer filter? I was wondering if any one with experience with both sensors can give some insight on this.

I've heard there can be problems with alignment of the different signals using a 3CCD.

However, when using the Bayer filter, you need to convert, and you have problems of lower resolution, introducing colors that are not there in reality, etc.

Is one of the two best usable? Or is it application specific?

EDIT:
Maybe I can add some extra questions to make the question more clear. For the record, I'm not asking this question for a specific project, I'm just interested in the differences between 3CCD or Bayer in general. I know the basic (hardware) differences between 3CCD and Bayer, but I'm looking for general details about the differences.

These type of questions come to my mind:

In an application that is based on color will it always add significant color information when using 3CCD instead of Bayer?
Are there other solutions instead of 3CCD to overcome the problems with Bayer? (thanks Ktuncer)
What exactly are the significant differences between 3CCD and Bayer (not application specific, and if application specific, in what kind of applications could it be significant)?

Generally Bayer mask (also known as CFA = color filter array) imagers are the most popular.

3Chip cameras cost considerably more - you are not only paying for 3x as many chips but the alignment in manufacturing is expensive and the range of cameras is much more limited. The largest are HD (1920x1080) although a lot of cameras actually have much smaller chips and rely on a deliberate half-pixel misalignment to interpolate each colour into a higher resolution - as with a CFA. Since most are aimed at the TV market they are usually video output (so you need a grabber card) and are a fixed TV frame rate

There are potentially aliasing issues with CFA cameras if you have a very small point image which may hit only one of the colour cells, this can be fixed with an anti-alliasing filter which blurs the image slightly, or by software.

One area in which 3chip wins is in low light where you need high colour quality. A 3chip camera detects essentially all the light in the scene in each colour. A CFA will detect only 25% of the Red and Blue and 50% of the green since the other photons are blocked by filters.


Why Does the DFT Assume the Transformed Signal Is Periodic? In many signal processing books, it is claimed that the DFT assumes the transformed signal to be periodic (and that this is the reason why spectral leakage for example may occur).

Now, if you look at the definition of the DFT, there is simply no that kind of assumption. However, in the Wikipedia article about the discrete-time Fourier transform (DTFT), it is stated that

  When the input data sequence $x[n]$ is $N$-periodic, Eq.2 can be computationally reduced to a discrete Fourier transform (DFT)

So, does this assumption stems from the DTFT?
Actually, when calculating the DFT, am I in fact calculating the DTFT with the assumption that the signal is periodic?

There are already some good answers, but I still feel like adding yet another explanation, because I consider this topic extremely important for the understanding of many aspects of digital signal processing.

First of all it is important to understand that the DFT does not 'assume' periodicity of the signal to be transformed. The DFT is simply applied to a finite signal of length $N$ and the corresponding DFT coefficients are defined by

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N},\quad k=0,1,\ldots,N-1\tag{1}$$

From (1) it is obvious that only samples of $x[n]$ in the interval $[0,N-1]$ are considered, so no periodicity is assumed. On the other hand, the coefficients $X[k]$ can be interpreted as Fourier coefficients of the periodic continuation of the signal $x[n]$. This can be seen from the inverse transform

$$x[n]=\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

which computes $x[n]$ correctly in the interval $[0,N-1]$, but it also computes its periodic continuation outside this interval because the right-hand side of (2) is periodic with period $N$. This property is inherent in the definition of the DFT, but it need not bother us because normally we're only interested in the interval $[0,N-1]$.

Considering the DTFT of $x[n]$

$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{3}$$

we can see by comparing (3) with (1), that if $x[n]$ is a finite sequence in the interval $[0,N-1]$, the DFT coefficients $X[k]$ are samples of the DTFT $X(\omega)$:

$$X[k]=X(2\pi k/N)\tag{4}$$

So one use of the DFT (but certainly not the only one) is to compute samples of the DTFT. But this only works if the signal to be analyzed is of finite length. Usually this finite length signal is constructed by windowing a longer signal. And it is this windowing which causes spectral leakage.

As a last remark, note that the DTFT of the periodic continuation $\tilde{x}[n]$ of the finite sequence $x[n]$ can be expressed in terms of the DFT coefficients of $x[n]$:

$$\tilde{x}[n]=\sum_{k=-\infty}^{\infty}x[n-kN]\tag{5}$$
$$\tilde{X}(\omega)=\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}X[k]\delta(\omega-2\pi k/N)\tag{6}$$

EDIT: The fact that $\tilde{x}[n]$ and $\tilde{X}(\omega)$ given above are a DTFT transform pair can be shown as follows. First note that the DTFT of a discrete time impulse comb is a Dirac comb:

$$\sum_{k=-\infty}^{\infty}\delta[n-kN]\Longleftrightarrow\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}\delta(\omega-2\pi k/N)\tag{7}$$

The sequence $\tilde{x}[n]$ can be written as the convolution of $x[n]$ with an impulse comb:

$$\tilde{x}[n]=x[n]\star \sum_{k=-\infty}^{\infty}\delta[n-kN]\tag{8}$$

Since convolution corresponds to multiplication in the DTFT domain, the DTFT $\tilde{X}(\omega)$ of $\tilde{x}[n]$ is given by the multiplication of $X(\omega)$ with a Dirac comb:

$$\begin{align}\tilde{X}(\omega)&=X(\omega)\cdot\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}\delta(\omega-2\pi k/N)\\&=\frac{2\pi}{N}\sum_{k=-\infty}^{\infty}X(2\pi k/N)\delta(\omega-2\pi k/N)\end{align}\tag{9}$$

Combining $(9)$ with $(4)$ establishes the result $(6)$.


Output of a system given it's transfer function and input (beginner) I have an exercise that gives me the transfer function of a system

$$H(s) = \frac{3s^2+27}{s^4+8s^3 + 16s^2} $$

and an input

$$x(t) = \frac13 cos(3t) $$

An ask's what is the output

I don't whant the answer but the steps to calculate it.

Where to start? I need to find the inverse or there's a shorcut?

For an LTI system, output $y(t)$ is given by
$$y(t) = h(t)\otimes x(t)$$
Where $x(t)$ is input and $h(t)$ is impulse response of the system. The operator $\otimes$ represents convolution.

Convolution operation is mapped into multiplication in Laplace domain. ie,
$$Y(s) = H(s)\times X(s)$$
Where, $Y(s)$, $H(s)$ and $X(s)$ are the Laplace transform of $y(t)$, $h(t)$ and $x(t)$ respectively.

You can use any one of the above equation to solve your problem.

Two ways:

Find $h(t)$ from $H(s)$ and use first equation. OR
Find $X(s)$ from $x(t)$, use second equation to find $Y(s)$ and then find $y(t)$ from it.

I prefer second method.


Acoustic Scenarios With Negative SNR For all the practicing engineers out there:

Where have you encountered a negative SNR?
How negative, in dB, was it or estimated to be?
Where would you expect to encounter a negative SNR?

Note: I am also counting interference signals as noise, so SNR can also include SIR.

In my own experience, I have seen negative SNRs in car communications. There is road noise, as well as multi-talker and music interference. I have also seen it in the output from cellular phones. Wind noise can sometimes pose a significant problem there, as well as the myriad of interferences in speakerphone mode.

I've heard rumors of covert acoustic room impulse response measurement being done during concerts. You would need the people in the room if you want to use the impulse response for artificial reverberation and you don't want the room to sound overly echoey like there's no-one in it. As received at the microphone, the signal would be a rather quiet noise-like test signal, and the music, the sounds from the people, and equipment noise would be the noise. Because the impulse response is in practice only a few seconds long, a repetitive test signal can be used. The recording would be accumulated into a circular buffer. Because noise has on average zero correlation with what has already been recorded in the buffer, its root mean square amplitude will grow on average as $\sqrt{N}$ for $N$ recording cycles, while the signal's root mean square amplitude will grow as $N$. At the end of the session, deconvolution is used to get from the test signal response to the impulse response. For a spectrally flat test signal such as a maximum length sequence (MLS), it suffices to convolve with the reverse of the test signal. For a buffer of $M$ sample points the deconvolution gain is $M$ for the signal and $\sqrt{M}$ for the noise. Combining the gains from repetition and deconvolution and assuming a sampling frequency of 44.1 kHz, we get a recording time $t$ dependent dB-scale gain difference of $20\log_{10}\left(\frac{44.1\text{ kHz }\times\ t}{\sqrt{44.1\text{ kHz }\times\ t}}\right)$ dB between the two:
Figure 1. Signal-to-noise ratio (SNR) gain for room impulse response measurement as a function of recording time.
Guessing some numbers for an "easy" scenario, if we start with the signal at a sound pressure level (SPL) of 30 dB and the noise at 60 dB SPL which is a starting SNR of -30 dB, in 1 hour of recording we get a SNR of 52 dB, growing very slowly after that, by 10 dB for every 10-times increase of the recording time. 52 dB is already useful for artificial reverberation.


Poisson Noise Removal from an Image What type of filters can be used to remove Poisson noise in an image? I tried applying Gaussian filter to few points in the spectrum.

In image denoising far more important then the noise distribution is the noise spatial correlation properties and the prior about the image.

Let's try building some cases and dealing with them.
The model is:

$ y = x + n $

Where $ x $ is the clan image, $ n $ is the Poisson Noise (With mean $ \lambda $) and $ y $ is the noisy image.

Noise Is Poisson Distributed and Spatially White

In this case each pixel $ {y}_{i} $ is distributed with Poisson with mean $ {x}_{i} + \lambda $.

Deriving the the estimator for $ {x}_{i} $ will yield it is $ {y}_{i} $. So no gain here what so ever.

Noise Is Poisson Distributed and Spatially White with Image Prior

If we assume some kind of a prior on the image (Piece Wise Smooth) we can do something.
The simplest idea is use the mean of a neighborhood we believe should have the same value (Think of capturing a wall with the same color all over it).

A very popular and effective way to enforce Piece Wise model is using the Total Variation prior on the Gradients of the image.
Then the MAP becomes something to work with (See Noise, Image Reconstruction with Noise (EE367/CS448I: Computational Imaging and Display, Class 10, Gordon Wetzstein, Stanford University)).

Noise Is Poisson Distributed and Spatially Correlated

The usual technique to deal with spatially correlated noise is working in a Multi Scale manner.

Noise Is Poisson Distributed and Spatially Correlated with Image Prior

Now combine all methods and voila, you have a real world denoiser.


Cross-Correlation Signal Delay Estimation Variance I am working on a project which intends to use Time Difference of Arrival (TDoA) for localization.

Firstly, my understanding is that a matched filter is the most common method for estimating signal-delay.

To get performance metrics for TDoA, I need the variance of the time-delay estimate. Thus, what I'm looking for is a way to calculate the variance in the signal-delay estimate using a matched filter as an estimator.

This seems like it would be a very common problem; however, I have consulted three books and searched online, and I have not found anything which addresses this issue.

One warning first: this is going to be a bit literature-heavy, because I think that's what you're more or less explicitly ask for. You want something that you can read, understand and cite for your publication.

  Firstly, my understanding is that a matched filter is the most common method for estimating signal-delay.

No. A matched filter can be used to effectively correlate a receive signal against a transmit signal shape. It's can but be a puzzle piece in a TDOA system.

  Thus, what I'm looking for is how to calculate the variance in the signal-delay estimate using a matched filter as an estimator.

Again, the matched filter is not the delay estimator; think of it as a signal preconditioner!

Now, what you probably mean is that you use matched filters for the TX pulse shape and the RX pulse shape filter, and then, based on the post-filter signal, do a delay estimate.

I'd like to cite a paper,[1] :

which doesn't actually deal with matched filtering directly, as the TDoA system isn't able to define the pulse shape the transmitter uses; instead, it's about cross-correlation of the received signals at different positions.

The interesting point is that the defining entity (aside from the obvious $\frac{c}{f_\text{Nyquist}}$ limit) is the SNR of the received signal -- that's pretty obvious, if you think about it: The "certainty" with which the correlation of RX and TX signal has a maximum at the "right" position is reduced with increased SNR.

In previous work of the author, she uses in [2] that with AWGN assumption and cross-correlation (i.e. matched filter with a known TX signal)

$$\begin{align*}
\sigma^2_{\hat\tau}&=\frac{3}{8\pi^2}\,\frac{1+2\, \text{SNR}}{\text{SNR}^2}\, \frac {1}{T_o B^3}\text,\,\text{with}\\
\sigma^2_{\hat\tau}&: \text{Variance of the timing estimate}\\
T_o&: \text{Duration of observation, and}\\
B&: \text{bandwidth of the system.}
\end{align*}$$

Interesting here, aside from the $B^{-3}$ relationship of variance and bandwidth, is the cited paper (if you find it in fulltext available online for free; I have IEEExplore access, but I don't know about you...) [9], which pretty much should give you plenty of background to work with.

[1]Noha El Gemayel, Holger Jäkel, Friedrich K. Jondral: Error Analysis of a Low Cost TDoA Sensor Network, IEEE/ION Position, Location and Navigation Symposium 2014, Monterey, CA, USA, May 5-8, 2014

[2]Noha El Gemayel, Sebastian Koslowski, Friedrich K. Jondral:
A low cost TDOA Localization System: Setup, Challenges and Results,
,
10th Workshop on Positioning, Navigation and Communication 2013 (WPNC'13), Dresden, Germany, March 20-21, 2013

[9] G. Carter, Coherence and time delay estimation, Proceedings of the
IEEE, vol. 75, no. 2, pp. 236 – 255, Feb. 1987


MMSE - How to minimize a complex error with respect to a set of real parameters Suppose there's a complex signal $X(k)$ (where $k \in \{0, 1, 2,...,N - 1\}$) corrupted by additive complex noise. Its estimate $\hat{X}(k)$ is a linear combination of a set of real parameters $A_r$ ($r \in \{1, 2, 3,..., R$})
$$\hat{X}(k) = \sum_{r=1}^R A_rZ_r(k)$$
where $Z_r(k)$ is complex (and known).

I wish to obtain the real values of $A_r$ for which the error is minimized. If I simply differentiate the mean squared error (MSE) $\frac{1}{N}\sum_{k = 0}^{N-1}(X(k) - \hat{X}(k))^2$ with respect to each $A_r$ and set the resulting derivatives equal to zero, the values of $A_r$ I'll obtain will be complex, so that's not a solution.

My question is: how do I obtain the optimal values of $A_r$ such that the MSE is minimized, under the constraint that each $A_r$ should be real?

A very common approach is to consider $X(k) \;\text{and} \;\hat{X}(k)$ as  elements vector space $C^R$, and consider the distance between the 2 vectors as a norm, which are real, positive or zero, and satisfy the triangle inequality.

So using vector-matrix notation with the vectors as column vectors:

$$ \text{error}^2 =(\mathbf{X}-\mathbf{\hat{X}})^H
(\mathbf{X}-\mathbf{\hat{X}})
$$
where $H$ is conjugate transpose.

Constraint is satisfied when imaginary part of $A_r$ is zero.

$$
\frac{1}{2}(A_r - A_r^*)=0
$$
which is appended to the objective as a lagrange multiplier(s)
$$
\text{error}^2+\sum \lambda_r \frac{1}{2}(A_r - A_r^*)
$$

These sorts of problems are easier when you use Brandwood derivatives,

D. H. Brandwood, "A complex gradient operator and its application in adaptive array theory," in Communications, Radar and Signal Processing, IEE Proceedings F, vol. 130, no. 1, pp. 11-16, February 1983.
doi: 10.1049/ip-f-1.1983.0003
Abstract: The problem of minimising a real scalar quantity (for example array output power, or mean square error) as a function of a complex vector (the set of weights) frequently arises in adaptive array theory. A complex gradient operator is defined in the paper for this purpose and its use justified. Three examples of its application to array theory problems are given.
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4645581&isnumber=4645575

This is the convention used in

Van Trees, Harry L. Optimum array processing: Part IV of detection, estimation and modulation theory. Vol. 1. New York, NY, USA: John Wiley & Sons, 2002.


Refining accelerometer noise using Kalman filter I have read that one of the uses of the Kaman filter is to refine the noisy sensor  measurements.
I am using an accelerometer to get the position by integration and want to use Kaman filter to refine the readings.
I know that the Kalman filter can be used in sensor fusion, but what are the steps when only one sensor is available and the physical model of the moving object is not available (the force or u vector is unpredictable).

but what are the steps when only one sensor is available and the physical model of the moving object is not available

Then it's not a Kalman filter.

A Kalman filter works because the system is observable.  In hand-wavy terms, you need to have redundant information about your system states, either because you have actual redundant inputs, or because you have an adequate* model of the system dynamics, and you're watching the system output over time.  If you can't do those two things somehow, then you don't have the prerequisites for a Kalman filter.

Taking GPS/IMU fusion as an example, your design process can be completely naive about the actual plant dynamics because you can model the plant as a rigid lump of stuff that takes rotation rates and accelerations as inputs, and gives geolocation as outputs.  Your model has states for position, velocity, and orientation, and because you're constantly correcting these against the GPS readings, you can get a pretty darned accurate estimate** of the states.

As another example -- and possibly as an inspiration for you -- I'm very slowly developing a Kalman filter for the motion of an object that's constrained to a 2-sphere (in less hoity-toity mathematics, that's the surface of a sphere in 3-dimensional space).  Because we have gravity, it'll know "down"; as soon as the object turns on the sphere it'll know its altitude -- but it'll never know the compass direction***.

So -- if you can somehow come up with redundant information, then you can possibly come up with a true Kalman filter.  This redundancy can be as tenuous as the one I'm using in my second example, where the fact that the object is always a fixed distance from a point tells me a lot.  It could be something like a button on an app that says "ok, Phone, I promise that I'm holding you still" -- I'm pretty sure that with a 6-degree of freedom IMU you could measure horizontal and vertical distances that way (i.e., you could make a tape measure app), but I doubt you could do more.

* What is adequate?  Enough to make the system observable.  Which is a circular definition within the context of this answer, because it would take too much space to explain.  It's a concept that you'll find about two weeks into a senior/graduate class in state-space control systems.  The bottom line for this answer is that it's a property of the system model -- and you don't have one of those.

** Assuming you have a rich enough motion profile -- in order to get good orientation information, the system needs to see varying acceleration, and the corresponding changes in position.  Such a system, if it was bolted to a rock, would be able to estimate position, velocity, and "down" very well, but unless it had gyros with errors significantly smaller than the Earth's rotation speed it wouldn't be able to tell you North and East.

*** Unless I use magnetometers, or get that really good gyro I mention above.


Equivalence of difference equation with system of difference equations I have this system and tried to write down its difference equation as follows

$$    y[n] = b_0v[n] + b_1v[n-1]\\
    v[n] = x[n] - a_1v[n-1]   $$

I was asked to prove that the system can be described by the following difference equation

$$    y[n] + a_1y[n − 1] = b_0x[n] + b_1x[n − 1]$$

Please help me how to prove it.

HINT:

Just express $y[n]$ and $x[n]$ in terms of $v[n-k]$, $k=0,1,\ldots$, and show the equality.


Impulse response measurement : practical questions To measure system impulse response I can do the following steps :

system input : it is excited with an excitation signal having a duration t_exc with a bandwidth bw_exc (for example sine linear sweep)
system output : it is measured during t_out with a sampling frequency f_s
transfer function H is obtained applying fft(out)/fft(in)
impulse response h is obtained applying ifft(H)

I have several questions about details of this method. I know that these questions really depends on system studied, but I would like to know the general guidelines :

How to determine the appropriate excitation signal duration?
Do you recommend excitation signal windowing?
How to determine appropriate output signal measurement duration?
before applying fft, how much zero padding do I need to add to excitation signal and output signal?
how to avoid zero division when excitation signal max frequency < f_s/2 ?

If this had to be the approach to measure the impulse response, here are my further suggestions specific to the qeustions:

How to determine the appropriate excitation signal duration?

The equivalent noise bandwidth of each FFT bin if not windowed is $1/T$ where $T$ is the duration of the signal in time prior to taking the FFT (note zero padding will not change this, but windowing will). So this will be the equivalent to the resolution bandwidth on network and spectrum analyzers, and will effect the resolution of the resulting transfer function in frequency. Resolution is not equivalent to number of samples, since we can zero pad in time which will result in more samples in frequency (over the same frequency duration) but this will just be a smoothing or interpolation of the same frequency domain result and will not help us further distinguish a closely spaced variation in the transfer function that may actually exist. To be able to measure such a tight variation in frequency, we need a longer time sample. Ultimately this is the expected length of the impulse response.

Do you recommend excitation signal windowing?

This depends on the excitation signal and the dynamic range of the resulting frequency transfer function. It would generally be a good idea but significantly complicates the processing; for example if you did end up using a swept sine as the excitation signal, the window will be significantly attenuating all frequencies represented at the start and end of the sweep, thus reducing the SNR and therefore quality of the results in the frequency domain transfer function for  those frequencies. Windowing decreases the frequency resolution of the main lobe of each DFT bin's equivalent filter at the benefit of significantly decreasing the sidelobes and thus improving overall dynamic range. Please see this post further detailing windowing considerations as this is what would all apply here:
How to calculate resolution of DFT with Hamming/Hann window?
Why would one use a Hann or Bartlett window?

How to determine appropriate output signal measurement duration?
Before applying fft, how much zero padding do I need to add to
excitation signal and output signal?

Both these questions affect the zero padding since we need the FFT's to be of the same duration prior to computing the frequency domain transfer function. The output signal will be the duration of the input signal plus the length of the impulse response. We already established that the input signal should be the duration of the input response, so in this case it would be zero padded out an additional amount. There is no harm in having a longer duration sequence, these are the minimums so if I were to do this it would likely be an iterative approach where I would evaluate both the sampling rate and signal duration from the frequency domain transfer function and time domain impulse response respectively: In the frequency domain, if the transfer function has not sufficiently attenuated by the time we reach half the sampling rate, aliasing will occur and the sampling rate should therefore be increased. In the time domain if the impulse response has not sufficiently attenuated by the time we reach the end of the duration of the resulting answer, time domain aliasing will also occur and the overall time duration of the signal should be increased.  To summarize with consideration to signal duration and zero padding: The signal duration must be longer in time than the impulse response duration. This signal must also be zero padded out further by the duration of the impulse response duration to allow for the output to fully settle. The output duration will match the input duration and no further zero padding will be required. The sampling rate must also be confirmed to be sufficiently high enough by reviewing the resulting frequency transfer function.

How to avoid zero division when excitation signal max frequency <
f_s/2 ?

Since you are creating the excitation signal you will know precisely if any frequency is 0, so this would be trivial to avoid. What I would be tempted to do personally is decide what dynamic range my solution needs to cover for the resulting frequency domain transfer function and then replace any frequency domain zeros with a very small number that I know is much smaller than that. However having an excitation signal that had lower energy in any given bin would be a poor choice for determining the frequency transfer function, which is required to be accurate to determine the impulse response. Granted in the design of this you would want to have an excitation signal with a bandwidth that sufficiently exceeds the bandwidth of the transfer function, in which case we could have its energy taper off at the higher frequencies that we know are well beyond the transfer function's bandwidth--and then it that case we would simply ignore (zero) the output beyond that frequency where the result is meaningless (and we wouldn't want it to be applied to the IFFT).
OTHER APPROACHES
There are a lot of drawbacks and challenges to determining an impulse response using this approach similar to the pitfalls in designing filters using the frequency sampling approach detailed here:
Difference between frequency sampling and windowing method
Alternate methods that would be more straightforward would be to slowly sweep a frequency tone in the time domain with a rate less than the resolution bandwidth of the measurement desired while comparing the amplitude and phase of the output and input directly in time (as done in a network analyzer). Alternatively but similarly to step through each frequency and compare the static results for each one.    Another approach would be to use a PRN sequence as the excitation tone or any other pseudo-random white noise source (random number generator) and use the Wiener-Hopf equations to determine the frequency domain transfer function which results in the least-mean-square solution for the transfer function. This is further detailed at this post with practical MATLAB code specific to doing this:
How determine the delay in my signal practically