Solving for impulse response h[n] given input-output pairs This is a homework problem that I've worked on and I want to confirm some of my reasoning.

I'm given two sets of input-output pairs of a particular system, $S$, that we know is linear time-invariant (LTI). Given those two pairs, find the impulse response h[n] of the system:

$$x_1[n] = [1, 0, 1] * S \longrightarrow y_1[n] = [0,1,0,2,0,1]$$

$$x_2[n] = [0, 1, 1] * S \longrightarrow y_2[n] = [0,0,1,1,1,1]$$

I found a linear combination of both inputs to yield the delta:

$$\delta[n] = x_1[n] - x_2[n-1] = [1, 0, 0]$$

Due to the system being LTI, the output to the delta is:

$$y[n] = h[n] =  y_1[n] - y_2[n-1] = [0,1,0,1,-1,0]$$

Since the input is 3 samples long and the output is 6 samples long, then my impulse response $h[n]$ has to be 4 samples long. So I discard the last two to finally yield:

$$h[n] = [0,1,0,1]$$

Convolving the input signals with the newly found $h[n]$ yields the given output signals. I am however uncomfortable with having to invoke the required length of $h[n]$ to yield its final form by deleting the last two samples.

Did I go wrong somewhere?

Without any generalizations, this particular problem could be solved with the following: Let the causal FIR impulse response of the LTI system be $h[n]$ of length $L$, which can be inferred from the first set that $L = 4$, hence let the impulse response be $$ h[n] = \{a,b,c,d\} $$

Then from the first set, it can be seen that
$$y[n] = h[n]+h[n-2] = \{a,b,c+a,d+b,c,d\} = \{0,1,0,2,0,1\}$$

from which you can deduce that $$a=0 ~,~ b=1 ~,~ c=0 ~,~ d=1 $$ which confirms that your result was right. For a longer impulse response you would need to solve a more involved set of (linear) equations...


3-tap FIR filter: simple expression for $H(e^{j\omega})$ using trigonometric identities We have a linear time-invariant system described by the input-output relation

$$y[n] = x[n] + 2x[n - 1] + x[n - 2]$$

Below is my approach to analyze this system.

The impulse response of this system $h[n]$ can be found if we input $δ[n]$ in $x[n]$'s position.

So
$$h[n] = δ[n] + 2δ[n-1] + δ[n-2]$$

From above, I can get an expression for $H(e^{j\omega})$ using the Fourier transform.
$$H(e^{j\omega}) = (1 + 2e^{-j\omega} + e^{-j2\omega})
= (1+e^{-j\omega})^2$$

My question is how to simplify this expression using trigonometric identities.

It would be pleasure if I can get some help.

HINT:

Write $H(e^{j\omega})=e^{-j\omega}G(e^{j\omega})$ and use $(e^{j\omega}+e^{-j\omega})=2\cos(\omega)$. $G(e^{j\omega})$ should turn out to be real-valued.


Digital filter simulating hydroacoustic signal distortion Preambule:
I'm designing a sound model for my small submarine game. Model is running on the server, and I want to present the client with a mono-channel wav-stream from his hydrophone (20kHz discretization should suffice, I target 20Hz-10kHz band). I want that signal to be relatively-realistic. Specifically, I need to account for the fact that water dissipates high frequencies much faster.

Sound model:
I plan to pre-generate (with help of some RNG and IFFT) raw propeller time-domain sound sample from it's intensity spectrum. Then I plan to amplitude-modulate it with shifted shaft rpm x propeller blade count sine wave to get that propeller beat. At this point i need to materialize that signal on the player's sensor. It will probably involve two filters:

Some filter (wich is the filter in question) wich will account for water dissipating signal energy.
Bandpass filter, wich corresponds to this particular sensor's sensitivity band.
Simple amplitude multiplication to account for range, sensor directivity etc.

According to R.J. Urik "Principles of Underwater Sound", on the distance of r meters from the target, band level (dB) in frequency bin f of such sound can be approximated using following formula:

[1]: ResultBL = SourceBL - 10 * log10(r * r) - r * F(f)
or without wave front expansion term (wich is trivial):
[2]: ResultBL = SourceBL - r * F(f)

where F(f) is some smooth monotonic function (big radical with multiple squared frequencies and constants) of frequency.

The questions (all about this one filter, or lack thereof):

Is it possible to design a time-domain digital filter for given range r, that implements\approximates the signal distortion [2]?
What type of filter would you recommend? Problem is not realtime, but calculations should be fast. I'm not limited to causal filters, so I'm sensing something like non-causal IIR?
What algorithms are applicable to this problem.
If such algorithm involves transfer function, how should I transform\approximate F(f) term in it?
Will it be fast? What factors affect it's speed?
Will such filter work good in the whole band (20Hz-10kHz)?
Will such filter work well on large intensity level variances, e.g. on both very weak and very strong signals?
Is such synthesis computationally-expensive?

Thank you.

I wouldn’t bother much with the precise shape of the filter because any real passive sonar will boost the higher frequencies at the receiver and any realistic source levels are going to around 70 years old and most likely a double screw, and aspect dependent. Cavitation is depth dependent. A single pole low pass filter with a corner around 10 Hz is probably ok.

Fidelity is one of those things that can become an obsession. There is ambient noise and self noise associated with hydro acoustic flow.  Higher frequencies exhibit hull shading.  No one band covers all the signal types of interest.

If you want to be obsessive, look at Mike Porter's Ocean Acoustics web site.

http://oalib.hlsresearch.com/

I'm partial to RAM PE and there at least was a broad band (decomposition of narrow band) demo for the MATLAB version.

There was a video game that came out around 20 years ago that collaborated with a company named SONOLYST.  They actually got an Oscar for sound effects for the Hunt for Red October movie.

I like Urick.  I actaully meet him a few times and took a week long class from him but you might find Ross a bit more helpful.

Mechanics of Underwater Noise, Ross, D. isbn={9781483160467},
  url={https://books.google.com/books?id=sdwgBQAAQBAJ, 2013, Elsevier Science


What Is Uncorrelated Noise In many applications such as estimation theory, when we need to estimate a parameter then we usually consider in presence of white gaussian noise of zero mean and some standard deviation. During Maximum likelihood estimation, we also use this assumption. So, my question is -

Do we consider noise to be uncorrelated or correlated in estimation?
What is the difference between correlated and uncorrelated noise and its significance
Why do we consider such a property of correlated or uncorrelated in estimation and when we say that measurement noise is gaussian.

1) It depends on the application, but in general, both. In some cases we can whiten the noise, say in the following situation:

$$Y = x+ \eta\text{ where }\eta \sim N(0, \Sigma)$$

where we want to estimate $x \in \mathcal{A}$ for some set $\mathcal{A}$ from $Y$. In this case, we multiply both sides by the inverse of the square root of $\Sigma$, and we have an estimation problem where $\Sigma^{-1/2} x \in \Sigma^{-1/2} \mathcal{A}$.

2) Two random variables are said to be correlated if their covariance is non-zero. In the cause of jointly Gaussian random variables, uncorrelated is equivalent to independent.

3) Uncorrelated things are often easier to deal with since things decouple (for example, when estimating in the MMSE sense, uncorrelated is a form of orthogonality (perpendicularity, geometrically)). Its the next best thing to independence (and in the most important case, the jointly Gaussian case, they coincide). As for when things are Gaussian, one may justify things via central limit theorem or just by saying the answer is usually easy to find, and several estimators coincide in the jointly Gaussian case (MAP, MMSE, MMAE) [In the jointly Gaussian MMSE case, the estimator is the linear MMSE estimator, which is quite easy to calculate and only relies on first and second order moments, which makes life a lot easier instead of relying on higher order moments].

MLE is just a paradigm for non-random parameter estimation (where as in Bayesian settings such as MAP, MMSE, MMAE, the quantity you are trying to estimate is modeled as random). None of these methods impose a constraint of signal+additive noise for parameter estimation.

You should try reading a book or two on the topic.


Grayscale morphology Hello I've been trying to understand the grayscale dilation operation using the following example, but my results are different.

$I = \begin{pmatrix}
\bf{5} & 2 & 4 \\
1 & -3 & 1
\end{pmatrix}$

is the image to be dilated, and the structuring element is

$S = \begin{pmatrix}
0 & x & \bf{0} \\
0 & 2 & -1
\end{pmatrix}$

I have marked the origin in the structuring element and the image (the bold $\bf5$ in the top left of $I$ and the bold $\bf0$ in the top right of $S$), and $x$ means that the pixel is not part of the structuring element.
The correct answer is supposed to be

$\begin{pmatrix}
5 & 2 & 5 & 2 & 4 \\
5 & 7 & 4 & 6 & 3 \\
1 & 3 & 1 & 3 & 0
\end{pmatrix}$

but my result is

$\begin{pmatrix}
4 & 7 & 5 & 6 & 4 \\
5 & 3 & 5 & 3 & 4 \\
1 & -3 & 1 & -3 & 1
\end{pmatrix}$

What I did was, put the center of the s.e. in each pixel, for each pixel in the neighborhood sum pixel to pixel their values and take the max of the neighborhood as the new value.

For example, when the s.e. is in pixel $I(0,1)$ then

$$
\max( I(0,0)+SE(0,-1), I(0,1)+SE(0,0), I(1,0)+SE(1,-1), I(1,1)+SE(1,0) ) =\\ \max(5+x=x, 2, 3, -4)
$$

Can someone explain to me what is the correct way of doing this?
In the dilation do we need to flip the structuring element? What about erosion?

According to this paper, the grayscale dilation of image $I$ by a non-flat structuring element $S$ is defined as follows:

$[I ⊕ S](x,y) = \displaystyle\max_{(s,t)∈S}\{I(x-s,y-t)+S(s,t)\}$

Since the origin of the structuring element is in the top right corner, we have that $s∈[0,1]$ and $t∈[-2,0]$ excluding the point $(s,t)=(0,-1)$. So using your example of element $I(0,1)$ we have the following:

$[I ⊕ S](0,1)=\max\{I(0,3)+S(0,-2),\ I(-1,3)+S(1,-2),\ I(-1,2)+S(1,-1),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ I(0,1)+S(0,0),\ I(-1,1)+S(1,0)\}$

$=\max\{0+0,0+0,0+2,2+0,0+(-1)\}=\max\{0,0,2,2,-1\}=2$

Which is the correct value for element $(0,1)$ in the dilated matrix. You can then simply apply the formula to the other $(x,y)$ value pairs of interest to obtain the final result.

The definition of erosion (and the other morphological operations) can also be found in that paper and applied in the same manner.


How to calculate SNR with sample noise and noisy signal only I want to calculate the SNR for a speech signal recorded within a noisy environment. I also sampled purely the noise within the environment.

The SNR is based on the ratio of power of clean signal and the power of noise signal. So my naive attempt to get SNR is:

$\frac{P_{signal+noise}}{P_{noise}} = \frac{P_{signal} + P_{noise}}{P_{noise}} = \frac{P_{signal}}{P_{noise}} + \frac{P_{noise}}{P_{noise}}= SNR + 1$

$SNR = \frac{P_{signal + noise}}{P_{noise}} - 1$

However I realized that:
$P_{signal+noise} \propto (A_{signal} + A_{noise})^2 = A_{signal}^2 + A_{noise}^2 + A_{signal}A_{noise}$
$P_{signal+noise} = P_{signal} + P_{noise} + avg(A_{signal}A_{noise}) \neq P_{signal} + P_{noise}$

I couldn't find anything about this in the Wikipedia article nor other threads with this particular question.

The key issue here is whether the signal and noise are uncorrelated.  Assuming two real random variables x and y are both zero mean, the power of the combined signal is

$E\{(x+y)^2\} = E\{x^2+2xy+y^2\} = E\{x^2\}+E\{y^2\}+2E\{xy\}$

If the two variables are uncorrelated, then $E\{(x+y)^2\}=E\{x^2\}+E\{y^2\}$  In other words, the power of the combined signal and noise is the sum of the signal power and the the noise power if the signal and noise are uncorrelated and zero-mean.  If this assumption is valid - as if often the case in practice - then the approach you outlined would work.  It would also mean that $P_{signal+noise}=A_{signal}^2+A_{noise}^2$.  If the signal and noise are correlated, then the approach would not work.

One thing you might consider is whether you are trying to measure the signal-to-noise ratio or the signal-to-noise-and-distortion ratio.  When the audio is present, the recorded signal will suffer distortion that will degrade the overall signal quality beyond what would be predicted simply via the noise power.  This would cause the SNR estimate obtained via the above approach to be optimistic in predicting the overall signal quality.


Bin sizes for non-uniform discrete Fourier transforms For a non-uniform discrete Fourier transforms, do the specified frequencies – i.e., $f_k$ in

– refer to the midpoint of the bin or the lower bound? I read the answer here, but that stated that the frequency bins for NUFFT are always evenly spaced which is not always the case for NUDFT (e.g., NUDFT-I and NUDFT-III – see here).

If they are the lower bound, does that mean $f_k$ includes frequencies $[f_k, f_{k+1})$?

The notion of a "frequency bin" is somewhat misleading. In general a sine wave of a single frequency will show up in ALL frequency "bins" of the Time discrete Fourier Transform. It will generate the highest value where $f_k$ is closest to the signal frequency and but the values at any other Fourier frequency are not zero, just smaller.

The only exception is if the signal frequency is an exact integer multiple of $f_k$ and not windowing is applied.


Power spectrum of downsampled image I am not sure if this is the right stackexchange community for this question, but here goes.

I am generating a Gaussian Random Field (GRF) of a pre-defined power spectrum. In the image below, the input power spectrum is shown (in each panel) as the solid, smooth red line.

As a sanity check, I get the power spectrum of the resulting realised GRF. In the image below, I show this power spectrum for the original, high resolution random field realization in the upper right panel in blue. As we can see, it follows the input spectrum closely, as required.

Now, for an application I want to downsample the original image (by a factor of 2; that is, every 2x2 block of pixels becomes one pixel in the new image). If I then take the power spectrum of this downsampled (and thus lower-resolution) image, it does not follow the input power spectrum anymore; this is shown in the lower-left panel.

My first instinct was that this might be a resolution effect. Therefore, I generated another realization of the GRF, but this time at a lower resolution (indeed the same resolution as the downsampled image). If it was a resolution effect, this power spectrum should display the same discrepancy. It does not, however: in the lower right panel you can see that this lower resoltuion spectrum also follows the input spectrum closely.

Which leads me to believe that the discrepancy is caused solely by the downsampling. This is puzzling to me however, since downsampling seems so straightforward (perhaps deceptively so).
Below, I also give the (python) code used to downsample the image.

So, rephrasing my question for clarity: Why is the power spectrum of the downsampled image not in qualitative agreement with the original power spectrum?

def downsample(self,im, fact=4):
    """
    Downsample an image. Default with 2x2 blocks, but this can be adjusted if desired
    Needed because first we generate at high res, then apply lensing,
    then downsample to Planck resolution
    """
    assert isinstance(fact, int), type(fact)
    sx, sy = im.shape
    X, Y = np.ogrid[0:sx, 0:sy]
    regions = sy/fact * (X/fact) + Y/fact
    res = ndimage.sum(im, labels=regions, index=np.arange(regions.max()+1))
    res.shape = (sx/fact, sy/fact)
    self.downsampled = res

I have solved the problem.

Since I was taking the sum of 4 pixels in downsampling, what I was doing in effect was convolving with a sum of delta functions (that is, one delta function for each pixel). In Fourier space this translates to multiplying with a cosine, which explains the increasingly strong (relative) discrepancy between input and observed power spectrum.

To mitigate this, Wikipedia helped me out. If we first implement a lowpass filter in Fourier space (filtering out high frequencies, to avoid aliasing later on), and then reduce the sampling rate (to 1 in 4 pixels), we get the desired result, and indeed: the observed power spectrum looks nice up until the cutoff:


How can I implement image widowing in my image in Java? My initial image is a medical image and my main objective is finding directionality. The original size of my image was 8X8 and I have rescaled the image to 32X32 to increase my accuracy. My input image is Fourier Transform Power Spectrum and my main objective is finding directionality. I want to window my image of size 32X32 as an image pre-processing approach. I want to implement the circularly symmetric window in Java. My main concern is that which window should I use and how can I implement? I tried to use the Gaussian window for my image. But, as I apply the Gaussian window, the whole meaning of my image changes. Am I doing something wrong?

I have implemented my Gaussian Window in Java as:

    double value = 0;
    double distance = 0;
    double x_center = (imageWidth - 1) * 0.5;
    double y_center = (imageHeight - 1) * 0.5;
    double cutoff = 2;
    double dc_level = 255;
    for (int i = 0; i < image.getWidth(); i++) {
        for (int j = 0; j < image.getHeight(); j++) {
            distance = Math.sqrt(Math.pow(Math.abs(i - x_center), 2)
                    + Math.pow(Math.abs(j - y_center), 2));
            value = dc_level * Math.exp(-1 * distance * distance)
                    / (1.442695 * cutoff * cutoff);
            powerspectrum[i][j] = (int) value;
        }
    }

Thank you so much for your help.

Consider the large image, and the region of interest (red):

If you will simply take the Fourier transform of the 8x8 region, the transform is assuming that the image is periodic like this:

The Fourier transform will see artificial directional features across the replicates of the 8x8 region. But if you pad the 8x8 image with zeros to 32x32, these artifacts disappear and the frequency resolution of the discrete Fourier transform is enhanced as a welcome side effect as well:

However, now a new problem is revealed: There are artificial directional features along the boundary lines of the 8x8 region, like in the upper half of the left boundary and in the left half of the upper boundary. A solution is to multiply the 8x8 image with a circularly symmetric window function:

The drawback is that the corners of the image are completely ignored. If you think it is useful, you could extend your region of interest so that you would not null the corners of that 8x8 region.

To test the performance of your dominant direction finding algorithm starting with the 8x8 input image, you can use test images. Start with a flat white 8x8 image. This should not give a dominant direction stronger than some level of approximation error. If it does, make sure that your window function has circular symmetry and that you have normalized fairly the calculation of all direction-specific integrals over the squared magnitude of the Fourier transform. Then try 8x8 images with directionality generated by setting the pixels to values $\cos(ax + by + c)$, where $\text{atan2}(a, b)$ defines the directionality, $x$ and $y$ are the pixel coordinates and $c$ should not affect the analysis result.

If your algorithm passes these tests it should be fine!


Textbook question regarding LTI System I was posed the following homework problem:

  2.10 The following input-output pairs have been observed during the operation of a time-invariant system:
  \begin{align}
x_1(n)&=\{\underset{\uparrow}{1}, 0, 2\}\overset{\mathcal T}{\longleftrightarrow} y_1(n)=\{\underset{\uparrow}{0}, 1, 2\}\\
x_2(n)&=\{\underset{\uparrow}{0}, 0, 3\}\overset{\mathcal T}{\longleftrightarrow} y_2(n)=\{\underset{\uparrow}{0}, 1, 0, 2\}\\
x_3(n)&=\{\underset{\uparrow}{0}, 0, 0,1\}\overset{\mathcal T}{\longleftrightarrow} y_3(n)=\{1,\underset{\uparrow}{2}, 1\}\\
\end{align}
  Can you draw any conclusions regarding the linearity of the system ? What is the impulse response of the system ?

I thought that I had a decent grasp on linearity and time-invariance, but I am having trouble relating them to a specific problem.  I know that linearity means that the system satisfies the superposition principle - and that the weighted sum of the input signals is equal to the weighted sum to each output signal.

How do I apply that to this problem?  I have the solution below, but I am unsure how they get there:

  2.10

  The system is nonlinear. This is evident form observation of the pairs
  $$x_3(n)\leftrightarrow y_3 (n)\text{ and }x_2(n)\leftrightarrow y_2(n).$$
  If the system were linear, $y_2(n)$ would be of the form
  $$y_2(n)=\{3, 6, 3\}$$
  because the system is time-invariant. However, this is not the case.

This looks like a very contrived question, but the idea is simple: if the system were linear, then, since $x_2[n]=3x_3[n+1]$, then $y_2[n]$ should be $3y_3[n+1]$. The impulse response can be obtained by time-shifting $x_3$ so that the $1$ is at $n=0$, and then time-shift $y_3$ the same amount.


Running an FFT filter on a large data set I have a non realtime application where I need to run a bandpass FFT filter on a data array of between 5k and 10k data points. Do I break it up into (say) 256 point chunks, run the FFT on that and simply move onto the next 256 items, or is that too naive?

You need to know the length of the impulse response of your bandpass filter, and use an FFT longer than 256 by that length (to prevent circular convolution boundary artifacts).  Combine the longer FFT results back into a sequence of length 256 windows by overlap-add or overlap-save.


convolution and frequncy response I have a sinusiodal signal x in time domain.I have done the convlotion as followings using matlab:
y=conv(x,x)  and y=conv(x,fliplr(x)).
Observing the result by Matlab, both of the above signal giving the same result in frequency domain. Why do you think is that?

P.S: FLIPLR(X) returns X with row preserved and columns flipped
    in the left/right direction.

X = 1 2 3     becomes  3 2 1
    4 5 6              6 5 4

Convolution in the time domain means you are multiplying their frequency transforms. This is one of the many properties of the Fourier Transform.

Now, when you flip a signal X[n], you can also write it down as x[n1], where n1 = -n. If you look at the time scaling property, try replacing 'a' by -1. You'll observe that the F(W) will now become F(-W).

Since F(-W) is exactly the same as F(W), your result does not change.


How to specify the input data bit width for fixed point implementation in windowed OFDM? I am implementing a fixed point model of Windowed OFDM and to quantize the input data, I used 16 bits(complex) to encode it. But I am not sure how do we specify the minimum number of bits that should be used to encode the input ensuring the output wont be "off" from the expected results.

I saw a similar question in Bit width for FFT but the person has said assuming "M" input bits. Is it an empirical data that you set by trial and error or is there some maths behind it.

Thank you

Depends on which processor platform you want to implement it. If it is an ASIC (Application Specific Integrated Circuit) chip, they come with ARM processor cores in most cases and real time computation require fixed-point implementation. There is also memory limitation in such processors due to smaller footprint of device. If you implementation is on processor which is 32 or 64-bit wide and has large amount of memory at its disposal (such as on Laptop, Desktop), you can use 32-bit fixed point notation (even though floating point would just work fine like IPP DSP Library from Intel). In your case, if your I and Q data each are 16-bit fixed point numbers, I don't think you have to worry too much. I am not sure if there is a specific advantage but almost all fixed point convention I have come across use 1 sign bit and rest fractional bit. In this case Q1.15 format. In that case you need to normalize all your symbols, channel co-efficients etc so that they are within +/-1.


Anti-Filter in Image Processing? Is there something like an anti-filter in image processing?

Say for instance, I am filtering an image using the following 13 tap symmetric filter:

{0, 0, 5, -6, -10, 37, 76, 37, -10, -6, 5, 0, 0} / 128

Each pixel is changed by this filtering process. My question is can we get back the original image by doing some mathematical operation on the filtered image.

Obviously such mathematical operations exists for trivial filters, like:

{1, 1} / 2

Can we generalize this to complex filters like the one I mentioned at the beginning?

Thanks, Anil.

What you're describing is call deconvolution. It's a concept frequently used for distorted images and for equalization in communications. As such you should be able to find a number of resources for your specific application. In general, the original image may not be recoverable exactly (it often isn't).

However, you can do a pretty good job depending on how severe the image has been undesirably filtered. In addition, many methods will attempt deconvolution by estimating the filter which it appears the signal has been affected by. If the assumption that the signal was convolved with a linear filter is a poor one, your performance may be severely inhibited.


How to use Robust PCA for image inpainting? I want to use Robust PCA for image inpainting, i.e. to estimate anomalous (bad) pixels. I assume that the bad pixels are known (mask).

I know how to calculate the Robust PCA using Alternating Direction Method of Multipliers but if somebody has already an implementation in Matlab I would not say no to sharing it. ;) But this is not my question.

How can I use Robust PCA for image inpainting, i.e. estimating the bad pixels? With Robust PCA I have X = L + S where X is my original image matrix for some training images, L is a low rank matrix and S is a sparse matrix.

As far as I can remember, one of the main interests of using R-PCA for image inpainting is that a pixel accurate defective pixel masks is not necessary.
Intuitively, the main idea is that the observed image has low rank, but that the defects (given by the mask);

disturb this low-rank property
form a sparse distribution.

Thus, image inpainting using R-PCA is given in this context by:

$X = L + S,$

where $L$ is the restored (inpainted) image and $S$ is the defective pixel mask.

This approach was implemented for example for background subtaction (motion detection task): the authors perform a low rank reconstruction of some shopping mall hallway. The moving pedestrians end up in the sparse component (=> mask) and the background is given by the low rank part. The backgournd does not have any "holes" in place of the pedestrians since a hole breaks the low rank constraint.

For more insight about this approach, see Yi Ma's low rank page.


How do you show the Discrete Time Fourier Transform of $x[n]=\cos(2\pi f_0n)$ is $ \frac{1}{2} \delta(f+f_0) + \frac{1}{2}\delta(f-f_0)$? How do you show the Discrete Time Fourier Transform of $x[n]=\cos(2\pi f_0n)$ is $ \frac{1}{2}  \delta(f+f_0) + \frac{1}{2}\delta(f-f_0)$?

Here is my thought process:

The definition of DTFT in my textbook is:
$$X(f) = \sum_{n=-\infty}^{\infty} x[n] e^{-j2\pi fn}\quad\text{for}\quad \frac{1}{2} \leq f \leq \frac{1}{2}$$

Therefore for this specific problem we have:
$$X(f) = \sum_{n=-\infty}^{\infty} \cos(2\pi f_0n) e^{-j2\pi fn}$$

I know from Euler's identity that
$$\frac { e^{i\theta} + e^{-i\theta}}{2} = \cos\theta$$

So we can convert the cosine above into the complex exponential like so:
\begin{align}
X(f) &= \sum_{n=-\infty}^{\infty} \frac { e^{j2\pi f_0n} + e^{-j2\pi f_0n}}{2} e^{-j2\pi fn}\\
&=\sum_{n=-\infty}^{\infty}   e^{-j2\pi n(f-f_0)} + e^{-j2\pi n(f+f_0)}
\end{align}

This looks very close to the $\delta(f+f_0) + \frac{1}{2}\delta(f-f_0)$ I'm looking for - but I'm not sure how to get there.

The sequence $ ~x[n] = \cos(\omega_0 n) ~$ , $-\infty < n < \infty$, is neither absolutely nor square summable, therefore its DTFT formally does not exist; i.e., the DTFT sum does not converge to a finite number, but diverges to infinity.
Because of the extreme importance of that signal in the context of both the theory and the practice of signal processing, we would like to express its DTFT in an adequate and consistent manner under all operations of signal processing. One such possible representation makes use of the Dirac impulse functions with the following observation:
Let $x[n] = 1$, for all $n$, be a discrete-time (unstable) signal whose DTFT  $X(\omega)$ given by the sum: $$ X(\omega) = \sum_{n=-\infty}^{\infty} 1 e^{-j \omega n} \tag{1} $$
This sum, resembles a CTFS (continous time Fourier series) synthesis relation, remembering the fact that the DTFT $X(\omega)$ is a periodic function of continuous -frequency $\omega$. Then the above sum in Eq.1 is re-interpreted as a CTFS synthesis as:
$$ X(\omega) = \sum_{n=-\infty}^{\infty} a_n e^{-j \frac{2\pi}{W} n \omega } \tag{2} $$ where  $\frac{2\pi}{W}$ is the fundamental frequency with the fundamental period being $W = 2\pi$, and the CTFS coefficients are all $a_n = 1$. From CTFS pairs, we find the signal $X(\omega)$ whose CTFS coefficents are all $1$, as the Dirac impulse train :
$$ X(\omega) = 2\pi \sum_{n=-\infty}^{\infty} \delta(\omega - 2\pi n) \tag{3}$$
and using the CTFS analysis equation, we have :
$$
\begin{align}
a_n &= \frac{1}{2\pi} \int_{-\pi}^{\pi} \left( 2\pi \sum_{n=-\infty}^{\infty} \delta(\omega - 2\pi n) \right) e^{-j \frac{2 \pi}{2\pi} \omega n } d\omega \\\\  &= \frac{1}{2\pi} \int_{-\pi}^{\pi} 2\pi \delta(\omega) e^{-j \omega n } d\omega = 1 \tag{4}
\end{align}
$$
The integral in Eq.4 is identical to an inverse DTFT of the impulse function $2 \pi \delta(\omega)$, hence we deduce that the DTFT $X(\omega)$ of $x[n]=1$ is
$$ x[n] = 1 \implies X(\omega) = 2\pi \delta(\omega) \tag{5}$$
Having found the DTFT of the sequence $x[n]=1$, we apply an estabslihed property of the DTFT :
$$ x[n] \leftrightarrow X(\omega) \implies x[n]e^{j\omega_0 n} \leftrightarrow X(\omega - \omega_0) \tag{6}$$
$$ 1 \leftrightarrow 2 \pi \delta(\omega) \implies 1 \cdot e^{j\omega_0 n} \leftrightarrow 2\pi \delta(\omega - \omega_0) \tag{7} \\\\$$
to the seqeunce $x[n] = \cos(\omega_0 n)$ after decomposing it into complex exponentials:
$$
\begin{align}
x[n] &= \cos(\omega_0 n) = 0.5 \left( e^{j \omega_0 n} + e^{-j \omega_0 n} \right)  \\ \\
X(\omega) &= \pi \delta( \omega - \omega_0) + \pi  \delta(\omega + \omega_0)
\tag{8} \\ \\
\end{align}
$$
Eq.8 gives us the requested DTFT of the sequence $x[n] = \cos(\omega_0 n)$.


Calculating the PDF of a waveform from its samples A while ago I was trying different ways to draw digital waveforms, and one of the things I tried was, instead of the standard silhouette of the amplitude envelope, to display it more like an oscilloscope.  This is what a sine and square wave look like on a scope:

The naïve way to do this is:

Divide up the audio file into one chunk per horizontal pixel in the output image
Calculate the histogram of sample amplitudes for each chunk
Plot the histogram by brightness as a column of pixels

It produces something like this:

This works fine if there are a lot of samples per chunk and the signal's frequency is unrelated to the sampling frequency, but not otherwise.   If the signal frequency is an exact submultiple of the sampling frequency, for instance, the samples will always occur at exactly the same amplitudes in each cycle and the histogram will just be a few points, even though the actual reconstructed signal exists between these points.  This sine pulse should be as smooth as the above left, but it isn't because it's exactly 1 kHz and the samples always occur around the same points:

I tried upsampling to increase the number of points, but it doesn't solve the issue, just helps smooth things out in some cases.

So what I'd really like is a way to calculate the true PDF (probability vs amplitude) of the continuous reconstructed signal from its digital samples (amplitude vs time).  I don't know what algorithm to use for this.  In general, the PDF of a function is the derivative of its inverse function.

PDF of sin(x): $\frac{d}{dx} \arcsin x = \frac{1}{\sqrt{1-x^2}}$

But I don't know how to calculate this for waves where the inverse is a multi-valued function, or how to do it fast.  Break it up into branches and calculate the inverse of each, take the derivatives, and sum them all together?  But that's pretty complicated and there's probably a simpler way.

This "PDF of interpolated data" is also applicable to an attempt I made to do kernel density estimation of a GPS track.  It should have been ring shaped, but because it was only looking at the samples and not considering the interpolated points between the samples, the KDE looked more like a hump than a ring.  If the samples are all we know, then this is the best we can do.  But the samples are not all we know.  We also know that there's a path between the samples.  For GPS, there's no perfect Nyquist reconstruction like there is for bandlimited audio, but the basic idea still applies, with some guesswork in the interpolation function.

What I'd go with is essentially Jason R's "random resampler", which in turn is a presampled-signal based implementation of yoda's stochastic sampling.

I've used simple cubic interpolation to one random point between each two samples. For a primitive synth sound (decaying from a saturated non-bandlimited square-like signal +even harmonics to a sine) it looks like this:

Let's compare it to a higher-sampled version,

and the weird one with the same samplerate but no interpolation.

Notable artifact of this method is the overshoot in the square-like domain, but this is actually what the PDF of the sinc-filtered signal (as I said, my signal is not bandlimited) would also look like and represents the perceived loudness much better than the peaks, if this were an audio signal.

Code (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -> [Double] -> [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)

rand list is a list of random variables in range [0,1].


Computation Costs of Filters Is there any general knowledge regarding the costs of computation time (computation complexity) of filters? Like what types are "generally" less/more costly than others.

it's in textbooks.  you'll find that, for a sharp filter of given sharpness, that IIR filters will cost less than FIR.

given it's IIR, you'll find that the Direct Form I and Direct Form II will have less computational cost than the State-Variable Filter, Lattice, Normalized Ladder, Gold-Rader, or Zölzer, or Harris-Brooking forms (and there are a zillion other forms, i s'pose, but i can't see them being cheaper to implement than the Direct Forms).  and then you can implement many forms (like DF1, DF2, or Lattice) as one-big direct filter section or as a cascade of second-order sections.  but there are numerical tradeoffs doing it the cheap (from a computation time POV) way.

given it's FIR, there are at least 3 different ways of implementing an FIR, all with different cost issues.  if the FIR filter is totally general (but still causal), then the choice is between the simple direct method (the "transversal" FIR form) and what we like to call Fast Convolution (implementing a long FIR using the FFT).  for some FIRs of a greatly restricted class (that is FIRs that look like IIRs), there is a method that is called Truncated IIR (TIIR) that is about as cheap as an IIR, but it's mathematically FIR.

good question, no short answer that is also comprehensive.  need to get O&S and read it.


Confused about 1D log Gabor I am doing offline signature verification for my project. While researching online I found that tha 1D log Gabor method can be used to find Hamming distance of an image. I got stuck in understanding 1D log Gabor. Wikipedia have some information. My question is how can we find the value of $f$, $f_0$ and $\sigma$ of an image given in Wikipedia 1D formula ?

Here is the link where I got idea about log Gabor can be one method for verification.

You're misunderstanding:

Wikipedia describes $G(f)$, the frequency response of the Log-Gabor filter. So you don't "find" $f$, it's your free variable.
$f_0$ and $\sigma$ are parameters of your filter, which you chose to maximize the "distinguishability" of "different" pictures. You need to come up with what "different" and "similar" means for you, and then you can optimize these parameters to make the Gabor transform of picture reflect these definitions. Notice that this is totally dependent on the image content and application you're doing.


FFT on circularly symmetric function I have an circularly symmetric image in Fourier space:

where the 0-frequency is in the center. Now, I want to take the FFT of this image. Having numpy's FFT conventions in mind (i.e. 0 frequency at (0,0)), I first apply the the inverse-fftshift, and then apply the inverse-fourier transform. As follows:

image = np.fft.ifft2(np.fft.ifftshift(FT_image))

This results in the following real space image:

But I can't really wrap my head around this result. It looks as if another fftshift would be in order, but in real space this shouldn't be necessary.
I would expect a radially symmetric transform, centered in the center of the image.

Is this expectation wrong, or is my computation wrong? And either way: can you explain why?

Your Fourier space image is purely real, so the inverse Fourier is going to be radially symmetric about the origin. Because the origin is at (0,0), it takes a shift to put the inverse image in the center of your matrix.


What is the Implementation of a Low, Mid or High Tone Dodge or Burn in GIMP? In The GIMP, and also in photoshop, you have the DODGE/BURN tool to enhance colors in shaded or highlighted regions, which would otherwise be blackish or whiteish. I found empirically that this filter is a good preprocessing step for a computer vision task I am implementing.

I found here that there are 6 formulas which simply map RGB values to other RGB values independently of channel or pixel location. According to the link the transformation is y=x+0.25*sin(x) for midtone dodging.

However, I could not verify this with the gimp source code, an excerpt from gimpoperationdodgemode.c:

gboolean
gimp_operation_dodge_mode_process_pixels (gfloat              *in,
                                      gfloat              *layer,
                                      gfloat              *mask,
                                      gfloat              *out,
                                      gfloat               opacity,
                                      glong                samples,
                                      const GeglRectangle *roi,
                                      gint                 level)
{
  const gboolean has_mask = mask != NULL;

  while (samples--)
    {
      gfloat comp_alpha, new_alpha;

      comp_alpha = MIN (in[ALPHA], layer[ALPHA]) * opacity;
      if (has_mask)
    comp_alpha *= *mask;

      new_alpha = in[ALPHA] + (1.0 - in[ALPHA]) * comp_alpha;

      if (comp_alpha && new_alpha)
    {
      gint   b;
      gfloat ratio = comp_alpha / new_alpha;

      for (b = RED; b < ALPHA; b++)
        {
          gfloat comp = in[b] / (1.0 - layer[b]);
          comp = MIN (comp, 1.0);

          out[b] = comp * ratio + in[b] * (1.0 - ratio);
        }
    }
      else
    {
      gint b;

      for (b = RED; b < ALPHA; b++)
        {
          out[b] = in[b];
        }
    }

      out[ALPHA] = in[ALPHA];

      in    += 4;
      layer += 4;
      out   += 4;

      if (has_mask)
    mask++;
    }

  return TRUE;
}

The essence of the code would boil down to these two lines. The other code is for handling transparency of the brush and iterating over pixels and color channels.

gfloat comp = in[b] / (1.0 - layer[b]);
comp = MIN (comp, 1.0);

Suggesting a general formula for all tones and both burn and dodge of  y = x / (1- f(x)). I could not trace the value of layer[b] back to anything sensible, and hence do not know what f(x) is supposed to be.

Moreover, the Wikipedia page on dodge and burn, and general pages about tone mapping, only mention the principle, but not any formulas.

Are there different definitions of dodge and burn filters? What is the definition, in the form of a formula y=g(x), of the gimp tool?

It basically using what Photographers calls "Shadow Mask" / "Highlight Mask".

If the image is in the range [0, 1] "Highlight Mask" is actually the value of the pixel.
"Shadow Mask" is the inverse, namely 1 - Highlight Mask.

Use this mask as linear interpolator of the result with the originals.


Signals and Systems - LTI - Transforms - Impulse Response I have $x(t)-> LTI -> y(t)$ where $LTI=h(t),H(jw)$.

and

$H(jw)=ab/((a+jw)(b+jw))$ where a and B are real numbers.

I am wanting to find the impulse response $h(t)$ as well as the input/output differential equation to the system.

If my understanding is correct I would want to take the Fourier transform of the frequency response/input and I am not sure where to go from there.

$$H(jw) = \frac{output}{input} = \frac{Y(jw)}{X(jw)}$$

$$\frac{Y(s)}{X(s)} = \frac{ab}{(a+s)(b+s)}$$

$$s=jw$$

Then you can use some of the notions from this example to solve for your impulse response / differential equation. Hope that helps!


Voss algorithm and pink noise the voss' algorithm generated a 1/f series. What does the series mean?
For example, we use three dice and the result sequence length is 8. Let's set the sequence is 8,15,6,9,14,10,9,7. Then what does this sequence mean? does it mean the sequence of frequency or something else?

The sequence produced by the Voss algorithm is a sequence of random numbers whose power spectral density decays with frequency as $1/f$. In other words, if you find the magnitude spectrum of the sequence (in other words, the magnitude of the DFT of the sequence), you'll see that it decays roughtly as $1/f$.

As to the meaning of the sequence: it doesn't have any intrisic meaning. One way to think about it is as noise with a particular frequency spectrum. In other words, if you normalize the sequence and feed it to your sound card, you'll hear a certain kind of noise. This noise can be added to other signals, such as audio, to give it specific properties.

If all you need is $1/$ noise, there are other methods to generate it; see Pink ($1/f$) pseudo-random noise generation.


Prove that convolution with a separable filter is equivalent to convolution on each 1D filter I would like to prove that convolution of an image $I \in \mathcal{M}_{m_1 \times n_1}$ with respect to a separable 2D filter $F$, (i.e., $F = F_1 F_2$, where $F \in \mathcal{M}_{m_2\times n_2}(\mathbb{R})$, $F_1 \in \mathcal{M}_{m_2\times1}$, and $F_2 \in \mathcal{M}_{1\times n_2}$) is equivalent to convolution with respect to $F_1$ and then $F_2$. If $A*B$ represents convolution of an image matrix $A$ under a filter matrix $B$, then we need to show that $I*F=(I*F_1)*F_2$.

We define convolution as such: $(I*F)_{(m,n)} = \sum_{i=-m_2/2}^{m_2/2}\sum_{j=-n_2/2}^{n_2/2}F[i,j]\cdot I[m-i,n-j]$, or $\sum_{i=-m_2/2}^{m_2/2}\sum_{j=-n_2/2}^{n_2/2} I[i,j]\cdot F[m-i,n-j]$.

Here's what I've tried starting with the left side (using the matrix indexing notation [m,n]):

$((I*F_1)*{F_2}){[m,n]} = \sum_{j=-n_2/2}^{n_2/2}(\sum_{i=-m_2/2}^{m_2/2}{I{[m-i,n-1]}} \cdot {{F_1}{[i,1]}}){[m-1,n-j]}{F_2}{[1,j]}$

Now, I notice that the terms ${F_1}{[i,1]}$ and $F_2[1,j]$ as placed in the summations can be combined to yield $F[i,j]$, however, we are indexing this entire term, so I don't believe we can combine these two terms as such. Additionally, the I don't see how to get terms involving $I$ to become $I[m-i,n-j]$.

Any advice would be much appreciated.

It requires a good deal of interpretation and a few clear sketches.I will try to provide a clean and easy approach, may be not very rigorous.

Let $x[n,m]$ be the original image (2D sequence) where $n$ is along the horizontal and $m$ is along the vertical axes and the origin being at the bottom left.

Consider the separable 2D filter sequence $h[n,m] = f[n]g[m]$ where $f[n]$ and $g[m]$ are the corresponding 1D filter sequences along the horizontal and vertical axes respectively.

The 2D output sequence of the convolution is denoted similarly as $$y[n,m] = x[n,m]\star h[n,m]$$ The procedure proceeds by writing the convolutiom sum:

$$ \begin{align}
y[n,m] &= \sum_k \sum_r h[n-k,m-r] x[k,r] \\
&= \sum_k  \sum_r f[n-k] g[m-r] x[k,r] \\
&= \sum_k  f[n-k] \left( \sum_r g[m-r] x[k,r] \right) \\
\end{align}
$$

At the this point, recognize the following: The summation inside the brackets indicate a 1D convolution between the filter $g[m]$ and the $k-th$ column of the original image: $x[k,m]$. Remember that during the summation the dummy index $k$ is fixed, whereas the other dummy index $r$ is ranging through the vertical samples of $x[k,r]$. We shall denote this intermediate sum as a new 2D sequence called $w[k,m]$, parametrized according to $k$ and $m$. That's the intermediate image defined for each value of $k$ and $m$ and formed by convolving each column-$k$ of $x[n,m]$ by the filter sequence $g[m]$. Then the equation follows as:

$$ \begin{align}
y[n,m]  &= \sum_k  f[n-k] \left( \sum_r g[m-r] x[k,r] \right) \\
&= \sum_k  f[n-k] w[k,m] \\
y[n,m] &= f[n] \star w[n,m] \\
\end{align}
$$

Where we interpret the last line as a 1D convolution between the rows of the intermediate image (for each fixed $m$) and the 1D sequence $f[n]$ along horizontal direction. Which finishes the proof.

Indeed, it's hard to visualize the process, as it's a two step approach. Nevertheless, relience on the assumed interpretations and verificaton by matlab coding would rather help.


number of cosine waves in DCT Why is it $N\cdot N$ a sufficient number of cosine waves equal to the size of the image or to a portion of it to describe it (or a portion of it) completely without losing information applying the DCT.

For example, in the JPEG encoding I divide the image into 8x8 blocks and I will have, for each block, a table of exactly 8x8 coefficients: one for each wave. But... why?

N orthogonal basis vectors can completely describe a point in N dimensional space.  Same for N * N orthogonal vectors and N * N samples describing an N*N dimensional space.

A DCT matrix is constructed from an orthonormal set of basis vectors.


Cross correlation vs FFT for finding phase between 2 signals I have been struggling with this issue for months now, trying to understand which one to pick: FFT or cross correlation.

I am given a sine wave (with phase, amplitude, frequency and time vector known), and I have to extract the phase shift between that sine wave and another signal. The new signal has the same frequency and the same time vector as the sine wave, however the amplitude and the phase are unknowns. What would you guys suggest? I am at a loss.

Thank you.

I don't see it as one or the other. You can use the FFT to more efficiently compute a cross-correlation (versus computing the cross-correlation in the time-domain).

There is lots of good information on this thread, about computing the "phase shift", from the offset of the peak in the cross-correlation function.


What does frequency weighting an acceleration signal mean? And How do I frequency weight the raw acceleration data collected? Disclaimer: I have very little knowledge in signal processing.
I'm trying to calculate a Whole Body Vibration(WBV) measure called Vibration Dose value(VDV).
$$VDV = \left(\int_{0}^{T}a(t)^4\,dt\right)^{1/4}$$
where $a(t)$ is a frequency weighted acceleration with weights $W_d$ for x, y acceleration components and $W_k$ for z acceleration component, over the frequency range 0.5 - 80 Hz.
I've read one literature saying since horizontal accelerations are more harmful, the weighting coefficients should be 1.4 for x, y and 1 for z axis - I do not know where this comes from.
Some literatures also show the following plot for frequency weightings used for analysis of acceleration signals.

Questions/Doubts:

I do not understand what the plot means and how to get weightings from it. But when I think about it, I feel like it has something to do with the acceleration data collected over 0.5 - 80 Hz and get the weightings from there. The thing is the weightings are not constant in the plot for the mentioned frequency region.

Also, how do I use the weightings with the raw-data collected to get frequency-weighted acceleration data?

Do I multiply the discrete raw-data with the the corresponding weights for each axis to get frequency weighted acceleration?

...I do not know where this comes from...

It looks like that's not a scientific article, so it doesn't have references at the end.  Read it carefully, and see if it references military or industry standards for vibration dosage.  Then read those standards and see if they refer back to any studies.
It'll be "fun", like black-water plumbing.  Think of it as the intellectual equivalent of cleaning a blocked sewer line.

I do not understand what the plot means and how to get weightings from it. But when I think about it, I feel like it has something to do with the acceleration data collected over 0.5 - 80 Hz and get the weightings from there. The thing is the weightings are not constant in the plot for the mentioned frequency region.

Probably you get the frequency spectrum of your accelerations in $x$, $y$ and $z$, then you multiply each of those by the individual weightings shown, then you calculate your $VDV$ from that.
I was going to say you could use Parseval's Theorem and just do the computation in the frequency domain -- but there's that $a^4$ term in there, so you may just need to weight it in the frequency domain, convert to time domain, and do the math.
But this smells strongly of being an established field, which means it has established methodologies; you should not invent your own method.   The established methodologies are often different from what an expert in signal processing would do, in detail if not in principle, and sometimes that different can be significant.  Even if your own special way is better in every regard, it'll still be "wrong" because it'll be different from how everyone else does it.
So find the source material via a literature search: either dig into the scientific literature (starting with the references in your second paper), or if you're trying to meet some regulation into the relevant industry or government standards.  Somewhere in all the muck the method will be spelled out.  If it's a regulation, it'll be written in techno-legal, and chances are you'll need to find an explanation of it someplace.  If you can't understand the actual defining document, then you can search on the title or designation of that document, i.e. "how do I implement SAE-234566-3?".

Also, how do I use the weightings with the raw-data collected to get frequency-weighted acceleration data?

Explained above -- convert your measurements into a spectrum, multiply said spectrum by the weightings in your graph.

Do I multiply the discrete raw-data with the the corresponding weights for each axis to get frequency weighted acceleration?

Yes.  Then, presumably, you take the weighted acceleration vector at each frequency point and you get it's Euclidian length, and it's that Euclidian length that you sum up.


How is 2D convolution calculated? I wish to implement the 2D convolution on an FPGA, so Ineed to understand how it is calculated in practice. The main difficulty that I found apparently 2 different ways showcases how to do it.
The first is to center the kernel to a position and calculate the sum of the overlapping products:
http://www.songho.ca/dsp/convolution/convolution2d_example.html
Making sense of indices in 2D convolution operations

And the other one that involves shifting the kernel
https://www.allaboutcircuits.com/technical-articles/two-dimensional-convolution-in-image-processing/

So what is the way?
Or perhaps they are equivalent in way that I don't recognise?

Technically, shifting the kernel above the still image, or shifting the image "below" the centered kernel are equivalent. This is because convolution of a kernel and an image is a commutative operation:

$$\sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} h(k_1,k_2)x(n_1-k_1,n_2-k_2) = \sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} x(k_1,k_2)h(n_1-k_1,n_2-k_2) $$

This however assumes infinite extend, etc. Whether it can easily be used to program an FPGA is a different topic.


How to eliminate noise that has the same frequency of my signal I have an audio signal , and I want to remove noise from it. I have applied a bandpass filter in order to remove frequencies that are not part of mi signal. But if there is background noise (some people talking) , it is not removed , because it has the same frequency range that my signal.
How can I filter my signal , to eliminate the noise that has the same frequency as my signal?

Whenever the desired signals and undesired signals are at the same
frequency band, you can't use the bandpass filter because it will remove both the desired and undesired signals.

In such cases statistical characterization for filter design is needed. Means you need to find out those statistical characteristics of background voice (noise) which differs your required signal from the noise. Examples include variance, expectation / mean or there can be others depending upon background noise, you might need to design adaptive filter if the noise characteristics changes with time.


FFT butterfly input index I am trying to determine a "simple" way to compute which inputs of a FFT need to "butterfly" together for its various stages. I am looking at a diagram like this:

In that diagram of a $16$ point FFT,

For stage 1: 0 butterflies with 8, 4 butterflies with 12, and so on.
For stage 2: 0 butterflies with 4, 8 butterflies with 12, and so on. Of course stage 2's inputs are stage 1's outputs.

My hope is that for each stage of the FFT, I can just have a simple counter going from $0$ to $N-1$ (FFT length $N$), and I can do very inexpensive bit manipulation to form the correct indices in the correct order.

Stage 1's pattern is fairly obvious: Have a 4 bit counter going from 0 to 15, and completely reverse the bits in the counter:

$0000 (0)$ with $0001\longleftrightarrow 1000(8)$
$0010\longleftrightarrow 0100(4)$ with $0011\longleftrightarrow 1100 (12)$
and so on.

I believe stage 2 is like this: bit reverse the counter, then swap the 2 leftmost bits.

Is there a generic formula for any stage of any FFT?

Assume $N=16$. so the counter $c$ chooses values between $0$ to $15$ and there are $s=1,2,\cdots,\log2(N)$ number of stages.
Based on your comment, the input-output relationship for stage $s$ can be given by:

Represent the counter in binary form in $\log_2 N$ bits (here 4 bits).
Circularly shift the counter $s-1$ times to the right.
Bit reverse the result.
Represent in decimal.

The following Matlab code can give the output in vector form as explained in the comment:

N = 16;
stage = 1              %Do it here for stage #1
c = 0:N-1;
c_bin = de2bi(c);
butterfly_idx = bi2de(fliplr(circshift(c_bin',stage-1)'))

Example for N=16 and stage=2`:

c_bin =

     0     0     0     0
     1     0     0     0
     0     1     0     0
     1     1     0     0
     0     0     1     0
     1     0     1     0
     0     1     1     0
     1     1     1     0
     0     0     0     1
     1     0     0     1
     0     1     0     1
     1     1     0     1
     0     0     1     1
     1     0     1     1
     0     1     1     1
     1     1     1     1

fliplr(circshift(c_bin',1)')

ans =

     0     0     0     0
     0     0     1     0
     0     1     0     0
     0     1     1     0
     1     0     0     0
     1     0     1     0
     1     1     0     0
     1     1     1     0
     0     0     0     1
     0     0     1     1
     0     1     0     1
     0     1     1     1
     1     0     0     1
     1     0     1     1
     1     1     0     1
     1     1     1     1

bi2de(fliplr(circshift(c_bin',1)'))

ans =

     0
     4
     2
     6
     1
     5
     3
     7
     8
    12
    10
    14
     9
    13
    11
    15

which is the expected result for the second stage.


Delaying output equals FIR filtering? I have a stream of audio blocks at input, each having 64 samples.

Sending a previous block of samples to output, i.e. delaying one block, means applying a unity gain FIR filter?

out[n] = 0 * in[n -  0] +
         0 * in[n -  1] +
         0 * in[n -  2] +
         ...
         1 * in[n - 64].

What's the frequency response of this filter?
Is there a use case for such design, or does this make no sense?

Yes this would properly delay the input by 64 samples but note that you can simply do:

$out[n]=in[n-64]$

Which was probably clear to you and you just wanted to show how it related to a 65 tap FIR with all taps = 0 except the last one.

The coefficients of the FIR filter represent the impulse response (strictly in discrete terms we call this the "unit sample response") of the filter. In this case the impulse response is a delayed impulse. The Fourier transform of a delayed (time shifted) impulse has a magnitude response in frequency that is all ones, and a linear phase response where the delay is equal to the phase slope according to:

$\text{delay} = d\phi /d\omega$  where $\phi$ is the phase in radians and $\omega$ is the radian frequency in radians/second.

We see this specifically from the Time Shift property of the Discrete-Time-Fourier-Transform:

$$y[n] = x[n-N]$$
$$Y(\omega)=e^{-j\omega N}X(\omega)$$

You can quickly plot the frequency response using Octave or Matlab as follows:

freqz([zeros(1,64) 1])

With the following result:

With a linear phase filter such as this, the signal is not distorted by the filter, only delayed. For more on that see this post: why is linear phase important


How to estimate the number of taps needed for subsequent channel estimation algorithms? I recently implemented an MMSE channel equalizer, and I have read about many other types of channel estimation/compensation algorithms out there. However, one thing that is never talked about it how to estimate the number of taps in your model of the channel you are trying to estimate. How is this accomplished?

For example, in the MMSE implementation I have, I know the training sequence, and thus, and use the training sequence to minimize my error. However, what if the channel impulse response is longer than the training sequence? I would never see its effects, and thus its effects on the training sequence will not be seen. How does one then estimate the number of taps needed for channel estimation?

Thanks!

The length of the impulse response is typically related to the frequency resolution of the channel transfer function. As a rule of thumb: the more detail there is in the frequency response, the longer the impulse response will be.

In practice there are a few things you can do: If you have full access to a similar, you can simply measure it with a very long impulse response measurement. Than you can truncate the impulse response and see what happens to the transfer function. The truncation will create errors and this way you can dial in the impulse response length to the point where the error is still tolerable.

You can also use physical knowledge about the channel. For example an audio amplifier has only a few electronic components, all of which are specifically designed to create a flat transfer function with little phase distortion. A handful of samples is fine for that. On the other hand  look at a loudspeaker in a room: the sound bounces around with multiple reflections until it finally dies off. In this case you would need many thousands of samples (not practical at all).

Many systems have bandpass or high pass characteristic: all acoustic systems are high pass because air can't transmit DC sound. Most communication systems are band pass since the information needs to stay away from the extreme edges of the band. In these case often the length of the impulse response is determined by the high pass roll off, i.e. the frequency and steepness of the high pass.


What happens, when an unlimited bandpass-signal gets mixed? In theoretical books there's always a bandpass-limited signal being mixed, that's been low-pass filtered beforehand.

But what happens, if the RF-signal from the antenna isn't band- or lowpass-filtered?

I know for example from decimation (lowpass-filtering and downsampling) that there would be aliasing, which isn't a form of mixing but nonetheless.

what happens, if the RF-signal from the antenna isn't band- or lowpass-filtered?

If the system is analog, lots of uncontrolled harmonics. If the system is digital, then yes, you will get aliasing.

Aliasing itself is not a form of mixing. But sampling is.

When an analog signal $x(t)$ is sampled at distinct times $\frac{1}{Fs}$ seconds apart, it is effectively mixed with a train of impulses at frequency $Fs$. Sample-and-hold turns these spikes into "steps". This is what is causing the "mirroring" of the spectrum.

In fact, by tuning your $Fs$ carefully, you can "shift" a signal into a baseband by undersampling it.

  I know for example from decimation (lowpass-filtering and downsampling) that there would be aliasing, which isn't a form of mixing but nonetheless.

Probably it is clear by now that it is. Take an analog signal $x(t)$, multiply it with a train of impulses $u_s(t)$ at some frequency $Fs$, now multiply the result by another train of impulses $u_d(t)$ that is equivalent to "taking out" every other pulse. If you don't first low pass $x(t)$ for the spectrum you are trying to "squeeze it in" post mixing by $u_s(t) \cdot u_d(t)$, you will get aliasing. It's the same concept applied twice.

Hope this helps.


Phase factors for an 32 point fft I have just started working on FFT and its related phase factors.

I have read for an N=2 point FFT, We have two phase factors, one on the positive axis and one on negative axis.

similarly for an 4 point FFT, we have 4 phase factors as theta=(pi/2), at

 0'  W 40=1
-90' W 41=-j
180' W 42=-1
 90' W 43= j

here we can see the values that are 180 degree apart are negative to each other.

Here I wanted to know how it works for an 32 point FFT.

As theta=(pi/16) do we have 16 phase factors on the positive axis and 16 phase factors on the negative axis with the interval of (pi/16) ?

so will it follow the pattern as I mentioned below?

 0'       W16 0 =1
-11.25'   W16 1 =0.98 -j0.195
-22.5'    W16 2 =0.92 -j-0.38
and so on
....
....
till
-168.75'  W16 15=
and the negative parts of the above values from **W16 16** to **W16 32**

Am I correct? please verify the above details and let me know how the phase factors works for an 32 point FFT

Thanks in advance

The base functions of the discrete time Fourier transform takes on values from the set $$\{z\in\mathbb{C}:z^N=1\}$$ where $N$ is the dimension of the vector space. The solutions can be parametrized to be in the form $$z_k = \exp(2 \pi i \frac{k}{N})$$ for $k=0..N-1$ because $$z_k^N = \exp(2 \pi i k \frac{N}{N}) = \exp(2\pi i k) = 1$$

In other words, the values the base functions components are chosen from are the $N$ non-identical $N$-th roots of 1, which are equidistantly distributed on the unit circle and include $1$ itself.


learnable segmentation or learnable edge detection Is there some learnable segmentation or learnable edge detection algorithm exist?

For example I mark N fotos by hands and than program "learns" how to do it best way and can do it on similar images.

there are some papers on this subject or source code?

Update:
I found some inforamtion supervised edge detection \ supervised segmentation
http://www.loni.ucla.edu/~ztu/Download.htm

but still looking for code.

There is a benchmark code and test and learning dataset to perform image segmentation with ground truth dataset at Berkeley website and the code in the same place.


T1 Mapping in Python? I'm trying to write a simple program to perform T1 mapping of a few MRI Images.
My understanding is that all I need is to plot the pixel intensity values by inversion time then fit a curve. The output of fitting the curve is what gives me the T1 values.
This was my attempt:
# Loop through every pixel in image
for row in range(rows):
    for col in range(cols):
        x = []
        y = []
        # Add all the inversion time and pixel intesities to list
        for image in dicom_images:
            x.append(image.InversionTime)
            y.append(image.pixel_array[row][col])
        # Fit curve
        t1 = np.polyfit(x, y, 2)

Assuming what I did above is right so far I am now stuck because I don't know how to take the three coefficients produced by np.polyfit and convert them into a single pixel value. I also don't know if there are additional steps I need to take.
I found some MATLAB examples online but they were really big and hard to follow. This is my first time doing T1 mapping so I'm trying to keep it simple.
Here are some plots of the data points:

You are fitting a polynomial to a logarithmic decay curve, which is probably wrong.  You probably want to fit a log curve to a log curve.  You can use the naïve method of np.polyfit(log(x), y, 1) which will probably work fine, or you can weight the points to make it more accurate, or you can use scipy.optimize.curve_fit.  All are described here: https://stackoverflow.com/q/3433486/125507
The coefficients returned by the polyfit are the coefficients of the polynomial.  So with deg=2, you're fitting a parabola to the decay curve, which again is probably wrong: p(x) = p[0] * x**2 + p[1] * x + p[2]  With the np.polyfit(log(x), y, 1) method you'll be fitting a straight line to a straight line, so there will only be a slope and an offset: p(x) = p[0] * x + p[1]  The slope would then be proportional to the decay rate of the log curve (somehow, which you can figure out for the type of value you need).
After getting this working, you can probably vectorize it for speed, fitting all pixels simultaneously instead of using for loops.  np.polyfit can handle this, for instance: "Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column."


simulation with SINR I have a DS-CDMA system and the chip width is fixed but every user have a different spreading code length. I want to have a simulation of signal to interference noise ratio (SINR) VS. bit error rate (BER). If I assumed that SINR is given like -25:0 dB. How can I calculate noise power (it is thermal noise)  and the interference power?

The SINR characteristics of a DSSS Multiuser Communication are very difficult to approximate theoretically.

For the Downlink at least, it is easy enough to characterize it with lesser parameters, but still difficult to obtain any particular derivation. It depends heavily upon the Spreading Sequences(their auto-correlations and cross-correlations), symbols of the interfering users(their average Energy) and so on. As a very decent approximation, $SINR$ is usually computed using random spreading sequences assumption. According to Andreas Goldsmith's textbook "Wireless Communications", for a $K$-user system with $N$ chips per symbol, the $SINR$(assuming Random signal transmissions for every user, again a probabilistic assumption) is given by $SINR = (\frac{N_{o}}{E_{s}}+\frac{K-1}{G})^{-1}$ where $G=N$ approximately.
For the Uplink it is even more a difficult task to analyze this. If all users' signals are received with the same power($AWGN$), it can be shown that for those Asynchronous users assuming Random Spreading codes with $N$ Chips per symbol, random start times and random carrier phases, $SINR = (\frac{K-1}{3N}+\frac{N_{0}}{E_{s}})^{-1}$.

For simple simulations, usually people assume that the $SINR$ can be calculated as, $SINR = \frac{E_{s}}{N_{0}+\sum_{i=1}^{K-1}\frac{E_{si}}{G}}$ where G is the processing gain which can be roughly taken to be as $K$. This factor of $G$ applies if the codes are perfectly orthogonal and its approximation to $N$ is usually justified assuming the properties of Spreading codes.

In all the above formulae, $N_{0}$ is the Noise Power which is equal to $2*\sigma^{2}$
for 2-Dimensional Signalling and $E_{s}$ is the average Received Energy of the signal. You can refer to Goldsmith if you to learn any further on the math behind it all. Now you can straight away plugin your numerical value on the LHS of the expressions that are stated above and you can go about calculating the Thermal Noise Power since you know the Received Powers $E_{s}$ for every user's signal. And all these formulae are for $AWGN$ conditions. The difficulty comes from the fact that Interferers cannot be approximated as Random Gaussian Noise.


Time invariance of a summation sequence The following sequence is given which is supposed to be time-variant:
$$y[n] = \sum_{k=n_0}^n x[k]$$
I'm having difficulties proving the time-variance or finding a counterexample for it being time-invariant.
My idea (which proves it being time-invariant?) is:
\begin{align}y_1[n] &= T\{x_1[n]\}\\x_2[n] &= x_1[n-n_0]\end{align}
check, whether    $T\{x_1[n]\} = y_1[n-n_0]$
$$y_2[k] = \sum_{k=n_0}^n x_2[k] = \sum_{k=n_0}^n x_1[k - n_0] = \sum_{k=n_0-n_0}^{n-n_0} x_1[k] = y_1[n-n_0]$$
Where is my mistake? Working around the summation is giving me trouble.
EDIT:
Following the advice in the comments for $n = n_0$ the left equation equals to 1 while the right equation equals to 0, making both equations unequal and therefore the sequence time-variant?
Second sequence
This sequence is supposed to be time-invariant.
\begin{align}y[n] &= \sum_{k=n-n_0}^{n+n_0} x[k]\\
y_2[k] &= \sum_{k=n-n_0}^{n+n_0} x_2[k] = \sum_{k=n-n_0}^{n+n_0} x_1[k - n_0] = \sum_{k=n_0-2n_0}^{n} x_1[k] = y_1[n-n_0]\end{align}
I followed the same steps but I'm still not sure about for the case $n = n_0$ since when added in left and right equation
$$\sum_{k=0}^{2n} = \sum_{k=-n}^{n}$$
They seem unequal to me?

Sequences that are "time-invariant" must be constants (not varying with time at all).
What the OP is actually asking about is a (discrete-time) system whose response $y$ to the input $x$ is defined to be
$$y[n] = \begin{cases}\sum_{k=n_0}^n x[k], & n \geq n_0,\\??? & n < n_0\end{cases}$$
where the $???$ is there because according to some traditionalists, if $n < n_0$, then that sum is an empty sum which has value $0$ but more liberal-minded might treat that sum as representing $$x[n] + x[n+1] + x[n+2] + \cdots + x[n_0-1] + x[n_0]$$ (with more serious-minded folks even negating that sum in analogy with the standard result that for $a<b$, $\int_b^a = -\int_a^b$). So, for the traditional meaning of things and $n_0$ denoting a (finite) integer, the system is not a time-invariant system. It is time-invariant (with unit pulse response $h[n]$ being the unit step function $u[n]$) if $n_0 = -\infty$ as noted in a comment by Robert Bristow-Johnson on the main question.


Calculating the fundamental frequency using Autocorrelation often gives half the expected value I'm currently writing a mobile app which needs to analyse musical notes and find the fundamental frequency to determine the pitch. To do this I'm reading in audio data, taking an FFT, taking the auto-correlation of that FFT and then finding the IFFT of that. Having done that I then find the maximum peak (which seems to always be a zero) and then find the next highest peak and determine the lag between them. Finally I calculate the period from this and find the fundamental frequency.

My problem is that this method is almost always returning a value which is half the expected fundamental frequency, only occasionally returning the correct value. Here is a graph of an autocorrelation that returned the correct fundamental frequency value:

When zoomed in on the relevant part we can see that there's a peak at around 370 which gives for the correct fundamental frequency of 117hz.

Here's a graph which returned the half the fundamental frequency value:

Zooming in on the relevant part we can see there's a trough at around 370 where we'd expect the peak to be:

Does anyone know why this happens and how I can resolve the problem/any other methods I can use?

...analyse musical notes and find the fundamental frequency to determine the pitch.

So, measure 440 Hz and output "A" (?), a tuner application.

  To do this I'm reading in audio data, taking an FFT, taking the auto-correlation of that FFT and then finding the IFFT of that.
  Having done that I then find the maximum peak (which seems to always be a zero) and then find the next highest peak and determine the lag between them. Finally I calculate the period from this and find the fundamental frequency.

This sounds like the superposition of two methods and two domains.

In the above process, once you obtain the FFT, you "transit" from the time domain to the frequency domain. In this domain, your signal is expressed as a sum of sinusoids at different frequencies.

Therefore, all you have to do (in this frequency domain) is locate the highest peak sinusoid.

This might be very easy if you are recording an electric instrument (bass, guitar, other) over wire or very challenging, if you are recording an instrument over the standard microphone found on a mobile phone.

BUT!, instead of using the FFT first, you could obtain the autocorrelation of the signal and then examine if that autocorrelation is periodic. The added benefit of this is robustness to noise.

Assuming that a dominant sinusoid exists in a recording, measuring the distance between the first and second peaks of the recording's autocorrelation would give an estimate of the period of the sinusoid. This technique applies to the time domain.

Because of the noise aspect, you might have to use a combination of the two. So, for example, record a signal, take its autocorrelation then apply FFT to the autocorrelation to discover its periodicity. In parallel to this, you might want to estimate spectral flatness as a measure of quality because it might be that the signal at the input does not contain any periodic component.

Finally, please note that due to the way discrete autocorrelation works, the autocorrelation of some signal $x(n)$ composed of $N$ samples, will have a length of $2 \times N - 1$ samples...which would need to be taken into account when trying to derive timings and frequencies from it.

Now, this is as far as what is mentioned in the question is concerned. In general, there are various techniques to perform pitch detection ranging from the very simple to rather complex ones. For a brief overview, please see this link, this link and this link, while, for something more specific, please see this link and this link

Hope this helps.


Comparison Between Guided Filter (Edge Preserving Filter) and Gaussian Filter I am looking for a good example that showing the different between Guided and Gaussian Filters. The example need to show the benefit of Guided filter, (for example: preserving edge...). Could you give me some example for that task? Thanks in advance

I tried some example, but it did not show benefit of Guided comparison with Gaussian

% example: edge-preserving smoothing
% figure 1 in our paper

close all;

I = double(imread('.\img_smoothing\cat.bmp')) / 255;
I = imnoise(I,'gaussian',0.1,0);
p = I;
r = 4; % try r=2, 4, or 8
eps = 0.2^2; % try eps=0.1^2, 0.2^2, 0.4^2

q = guidedfilter(I, p, r, eps);
std_Gb=1;
beta=0.1;
%% Initialization
Ng=ceil(3*std_Gb)+1; Gaussian = fspecial('gaussian',[Ng Ng],std_Gb);
imsm = conv2(I,Gaussian,'same');

[Gx,Gy] = gradient(q );
NormGrad = sqrt(Gx.^2 + Gy.^2);
Gb1 = 1./ (1 + 1* NormGrad.^2);

[Gx,Gy] = gradient(imsm);
NormGrad = sqrt(Gx.^2 + Gy.^2);
Gb2= 1./ (1 + 1* NormGrad.^2);

figure();
subplot(2,1,1);imshow([I, q,imsm],[]);
subplot(2,1,2);imshow([Gb1,Gb2],[]);

Create  a synthetic image with abrupt change form black to white (Better yest from Dark Gray to Bright Gray).
Add Gaussian Noise to it.
Filter it with Guided Filter and Gaussian Filter.
Look at the results.
Draw a 1D Line which crosses it and look at the Original, Noisy, Filters Guided, Filtered Gaussian.

Here are the results:

Noisy Image

Gaussian Filtered Image

Guided Filtered Image

As one could see, the Gaussian Filter smooth the image across high contrast edges while the Guided Filter keep the edges.

The code is located at GitHub - DSP StackExchange Q29041.


Why is my MATLAB's bode plot wildly off? I will give out all the details in case it is relevant. I have a MIMO state space system. I find its bode plot using MATLAB and separately using Mathematica. The plot from MATLAB is wildly off compared to the correct bode plot. However, the Mathematica plot is quite close to the correct one. Here's the interesting thing. The transfer function whose bode plot is being taken in both the softwares is calculated to be the same. Have I made some mistake or is it numerical errors contributing to the problem? Here is relevant part of my code:
A=[0,0,0,0,1,0,0,0;
    0,0,0,0,0,1,0,0;
    0,0,0,0,0,0,1,0;
    0,0,0,0,0,0,0,1;
    -297.3,163.5,0,0,0,0,0,0;
    162.9,-267.2,104.2,0,0,0,0,0;
    0,57.8,-74.2,16.4,0,0,0,0;
    0,0,16.4,-16.4,0,0,0,0];
B=[0,0,0,0,0;
    0,0,0,0,0;
    0,0,0,0,0;
    0,0,0,0,0;
    131.4,0.046,0,0,0;
    0,0,0.045,0,0;
    0,0,0,0.025,0;
    0,0,0,0,0.25];
S=[1,0,0,0,0,0,0,0;
    0,0,0,1,0,0,0,0];
D=zeros(2,5);
sys=ss(A,B,S,D)
systf=tf(sys);
s1a1=systf(1,2);
bode(s1a1);

[MATLAB plot which is way off the correct one][1][Mathematica Plot which is quite close to the correct one]2
EDIT: The transfer function being plotted here is:
$$\frac{0.046s^6-5.522\times 10^{-17}s^5+16.46 s^4 - 1.748\times 10^{-14} s^3 + 880.1 s^2 - 7.808\times 10^{-14} s + 7108}{s^8 - 4.373\times 10^{-15} s^7 + 655.1 s^6 - 2.258\times 10^{-12} s^5 + 9.887\times 10^{4} s^4 - 2.144\times 10^{-10} s^3 + 3.43\times 10^{6} s^2 - 4.353\times 10^{-9} s + 2.069\times 10^{7}}$$

I suspect that the reason for the difference is the number of samples used in the plot, since the extreme points are likely going to infinity (or very large numbers) such that if less points are used, those large values are simply missed in the computation. If the OP uses the same number of points in each, then the same result should be achieved- and for all the points shown if they were super-imposed on each plot it would be clear that for the points given, the results are the same.
My first suspicion before the OP added the plots was that it was due to the units of the axis. That is not the case but leaving this part of my answer below as it is another common reason for Bode plot differences:
The units of the frequency axis should be confirmed as they could be in units of:
Continuous domain units of frequency:
Hertz, or cycles/sec
radians/sec
Discrete domain units of (normalized) frequency:
cycles per sample (unique domain is $[-1,1)$)
radians/sample (unique domain is $[-\pi, \pi)$)


Finding the edge pixels on only one side of the object For my work on lung images, I get a binary image mostly resembling the pattern shown below. I need to find the pixels on the borders shown by yellow arrows. The edges on the other end are shown in pink color.

the edge function on MATLAB highlights all the edges irrespective of the side.

edge_BW_portionR = edge(BW_portionR)

I also used distance transform in which the lower right corner was used as a reference , and with respect to that, I put some conditions, but in vain.

How can this be solved?

I found a way to detect the right hand side edge with MATLAB:

1.- acquiring image

A=imread('001.jpg');
figure(1);imshow(A)

A1=A(:,:,1);
[sz1 sz2]=size(A1)

figure(2);[P,h1]=imcontour(A1,1);

this detects the contour, but function imcontour produces an additional small amount of points completely unrelated to the contour, and quite far away

hold all;
plot(P(1,:),P(2,:),'r*')  % check

2.- quik-fix for imcontour producing a few points far away.

Px=P(1,:);Py=P(2,:);
nx_out=find(Px>sz2)
P(:,nx_out)=[]
ny_out=find(Py>sz1)
P(:,ny_out)=[]
hold all;
figure(3);plot(P(1,:),P(2,:),'r*');axis ij  % check

3.- generating image with the complete contour only and 1 pixel thick contour only

Px=P(1,:);Py=P(2,:);
A2=zeros(sz1,sz2)
P=floor(P)
for k=1:1:numel(Px)
    A2(Py(k),Px(k))=255;
end
figure(4);imshow(A2)

4.- Selecting right hand side edge only

P2=[0;0]
for k=sz1:-1:1
    L1=A2(k,:)

    if sum(L1)>0
        nL1=find(L1>0)
        P2=[P2 [k;nL1(end)]];
    end
end
P2(:,1)=[]

figure(5);
imshow(A2);
hold all;

for k=1:1:length(P2)
    plot(P2(2,k),P2(1,k),'r*')  % check
end

The result can be improved refining the way pixels are selected, instead of

P2=[P2 [k;nL1(end)]];

using

P2=[P2 [k;nL1(2:end)]];

but this is assuming that the excluded left edge is a fairly straight up line, that may not be the case.

For more complicated shapes, to first find the single pixel thick contour
 have a look at my answer: https://uk.mathworks.com/matlabcentral/answers/350136-how-t-oget-straight-edge-in-binary-image?s_tid=srchtitle


What is the difference between PSOLA and TDHS time-scaling or pitch-shifting? What is the difference (if any) between PSOLA (Pitch Synchronous Overlap and Add) and TDHS (time domain harmonic scaling) time-pitch modification algorithms?

OK I did the two some time ago, TDHS in principle just apply time scale modification and to change the pitch do you need apply interpolation (resample) and it will shift the spectral envelope.

For TDHS is hard to find some paper that teach how its really works,  I learned the math and how it works in the Burazerovic Dzevdet paper:

$N_p$ is defined as the local Pitch Period

$\alpha$ is the time stretch factor

$N_c$ is the cross-fade length, here half the length of the window.

To compress the time $\alpha < 1$ the equation is:

$$ \begin{align}
N_c & = \operatorname{round}\left(N_p\frac{\alpha}{1-\alpha}\right) \\
& = \left\lfloor N_p\frac{\alpha}{1-\alpha} + \frac12 \right\rfloor \\
\end{align} $$

where "$\lfloor \cdot \rfloor $" is the floor() function returning the most positive integer not exceeding the argument.

The illustration from Dzevdet paper show how the algorithm works using a triangular window (n. b. 4 cycles become 3 and the period and fundamental frequency (pitch) of the waveform is not changed):

Now to expansion $\alpha > 1$ (note 3 cycles become 4 and the period and fundamental frequency (pitch) of the waveform is not changed):

$$ \begin{align}
N_c & = \operatorname{round}\left(N_p\frac{\alpha}{\alpha-1}\right) \\
& = \left\lfloor N_p\frac{\alpha}{\alpha-1} + \frac12 \right\rfloor \\
\end{align} $$

OK here my simple demo how TDHS works to compress:

alpha=0.5;
f=735;
Fs=44100;
signal= 0.9*sin(2*pi*f/Fs*(1:1000));
signal = signal';
period= 60;
out2=[];

Nc=round((period*alpha)/(1-alpha));
nsamples = length(signal);

ii=1;
out2=[];
while ( ii <= nsamples )

    if ii+Nc+period-1 > nsamples
        %OK I'm lost some samples in the end of out2 output signal, no problem just tests
        HI=000 %debug

    else

        frame = signal(ii:ii+Nc+period-1);
        frame1=frame(1:Nc) .* (1-linspace(0,1,Nc))';
        frame2=frame(period:Nc+period-1) .* linspace(0,1,Nc)';
        OUT = frame1+frame2;
        out2 = [out2' OUT']';

    end

    ii = ii + Nc+period-1

end

plot(out2)

Course this is just to show for you the concept, I used fixed integer period $N_p = 60 (44100/735)$

Now To PSOLA things will  different, Usually you need to use the Pitch information to mark the sign to the glottal instants you can change the Pitch and the time (no formants change here) .

Take a look in my answer about PSOLA here


The Sinc Function and the Partition of Unity Property For this question purposes a function, $ f \left( t \right) $ is said to hold the Partition of Unity if:

$$ \forall t, \ \sum_{n \in \mathbb{Z}} f \left ( t - n \right ) = 1 $$

Using the Poisson Summation Formula the Partition of Unity is equivalent of $ F \left( 2 \pi k \right ) = \delta \left ( k \right ) $ where $ F \left ( \omega \right ) $ is the Fourier Transform of $ f \left ( t \right ) $.

Looking at
$$ \frac{ \operatorname{sinc}^{p} \left ( pt \right ) }{\int_{-\infty}^{\infty} \operatorname{sinc}^{p} \left ( pt \right ) } $$
Does it holds the Partition of Unity for $ p \geq 2 $?

At first I thought it does, since the Fourier Transform equivalent property will hold.
Yet it seems I'm missing something.
Any ideas?

Thank You.

For any finite-energy signal $x(t)$, the function
$$y(t) = \sum_{n=-\infty}^\infty x(t-nT), -\infty < t < \infty$$
is a periodic function of $t$ with period $T$, that is, $y(t-T) = y(t), -\infty < t < \infty$. Note that $T$ is not necessarily the least period of $y(t)$.  Now, a periodic
function can be expressed as a Fourier series
$$y(t) = \sum_{m=-\infty}^\infty c_m \exp(j2\pi mt/T)$$ where
$$\begin{align}
c_m &= \frac{1}{T}\int_0^T y(t)\exp(-j2\pi mt/T)\,\mathrm dt\\
&= \frac{1}{T}\int_0^T \sum_{n=-\infty}^\infty x(t-nT)\exp(-j2\pi mt/T)\,\mathrm dt\\
&= \frac{1}{T}\sum_{n=-\infty}^\infty \int_0^T  x(t-nT)\exp(-j2\pi mt/T)\,\mathrm dt\\
&= \frac{1}{T}\sum_{n=-\infty}^\infty \int_{-nT}^{-(n-1)T}  x(s)\exp(-j2\pi ms/T)\,\mathrm ds\\
&=  \frac{1}{T}\int_{-\infty}^{\infty}  x(s)\exp(-j2\pi ms/T)\,\mathrm ds\\
&= \frac{1}{T}X\left(\frac{m}{T}\right)
\end{align}$$
where $X(f)$ is the Fourier transform of $x(t)$.
Now, for the case $T = 1$,
if $y(t) = 1$ for all $t$ as the OP desires, then it must be that $c_0 = 1$
and $c_m = 0$ for all $m\neq 0$, and therefore the Fourier transform of $x(t)$
must have the property that
$$X(m) = \begin{cases} 1, & m = 0,\\ 0, &m \neq 0,\end{cases}\tag{1}$$
that is, $X(f)$ is what is sometimes referred to as a Nyquist pulse (in the
frequency domain). We are familiar with Nyquist pulses in the time domain,
with $\operatorname{sinc}(t) = \frac{\sin(\pi t)}{\pi t}$ being the most
familiar one,  but there are many others, including as a special case
$[\operatorname{sinc}(t)]^k$. Applying duality,
$X(f) = [\operatorname{sinc}(f)]^k$ is a Nyquist
pulse in the frequency domain, and its inverse Fourier transform $x(t)$ is
a signal with the property that
$$\sum_{n=-\infty}^\infty x(t-nT) = 1, -\infty < t < \infty.$$
Note that $x(t)$ is the $k$-fold convolution of $\operatorname{rect}(t)$
with itself. For example, with $k = 2$,
$$x(t) = \begin{cases}1 - |t|, &-1 < t < 1,\\
0, & \text{otherwise.}\end{cases}$$
What about the other way around?  Well, $\operatorname{sinc}(t)$ has
Fourier transform $\operatorname{rect}(f)$ which has the Nyquist pulse
property $(1)$ in the frequency domain as does $[\operatorname{sinc}(t)]^2$
whose Fourier transform
$$X(f) = \begin{cases}1 - |f|, &-1 < f < 1,\\
0, & \text{otherwise,}\end{cases}$$
also satisfies $(1)$.  However, for $k > 2$,
the Fourier transform of $[\operatorname{sinc}(t)]^k$ does not
have the Nyquist pulse property $(1)$ because the three-fold (or more
generally $k$-fold convolution of $\operatorname{rect}(f)$ with itself
does not satisfy $(1)$: the $k$-fold convolution "bulges" outside
the interval $(-1,1)$.  See for, example, the
Irwin-Hall distribution
for details of the exact values of the convolution results.


Can we band-limit a PCM signal? Is it even relevant to relate a totally Analog term(Band limit) with digital one(PCM) Can we band-limit a PCM signal?

Yes -- it's just a matter of selecting the right pulse shape.

(By the way, it makes perfect sense to talk about the bandwidth of a PCM signal; it is, after all, a continuous-time physical signal with bandwidth, energy, power, etc, just like any other signal.)

Most introductory textbooks (irritatingly) assume that a PCM signal must use rectangular pulses. You don't have to, and actually you almost never do, except in the simplest applications.

The PCM signal is a train (or sequence) of pulses with varying amplitudes. This train can be written as $$s_{\text{PCM}}(t) = \sum_k a_k p(t-kT_p),$$ where $T_p$ is the pulse rate, and $a_k$ is the amplitude of each pulse.

In the simplest PCM signal, $a_k$ are either 0 or 1, and $p(t)$ is a rectangular pulse of duration $T_p$. The bandwidth of $s_{\text{PCM}}(t)$ is the same as the bandwidth of $p(t)$, which in this case is (mathematically) infinite.

By selecting a different pulse $p(t)$, you can limit the bandwidth of $s_{\text{PCM}}(t)$. A popular pulse shape is the "raised cosine", which has a bandwidth between $1/T_p$ and $1/2T_p$.

Note that the bandwidth can also be limited by "brute force": just filter the rectangular $s_{\text{PCM}}(t)$ with a low-pass filter of the desired cutoff frequency. This will distort the signal, but in simple applications (high SNR, no further distortion or attenuation) this usually causes no significant performance degradation.


Variance of a signal How to calculate the variance of noise samples modeled as follows:

$n_a(t)$ is a
Gaussian zero-mean white noise process with (two-sided) power spectral density $\frac{N_0}{2}$. $n_a(t)$ is passed through an anti-alisasing filter to get the output $n_f(t)$. How to calculate the variance of the noise samples $n[j]$ in terms of $N_0$ and $B$, where $n[j]$=$n_f(jT_s)$ and $T_s$ is the sampling period.

How to calculate the variance of the noise samples $n[j]$ in terms of $N_0$ and $B$, where $n[j]$=$n_f(jT_s)$ and $T_s$ is the sampling period?

Do you know how to calculate the variance of the process
$\{n_f(t) \colon -\infty < t < \infty\}$?
No? Hint: it is the area under the power spectral density curve of $\{n_f(t)\colon -\infty < t < \infty\}$ provided that you are not using radian frequency but are using Hertzian frequency to express the PSD.  Well, if you know the process variance, then be assured that each random variable $n_f(t)$ in the process has the same variance as the process variance. A more interesting question for you to ponder might be: "Are $n[k]$ and $n[\ell]$ independent random variables or correlated random variables?" which will bring in that nuisance parameter $T_s$ into play.


Finding squares in Image I need to find the squares in an image using OpenCV (no problem in matlab or any other, generally what i expect are some ideas).

Consider the test image below :

I need to find those coloured squares in above image accurately (not the white long strips).

What I have done :

I applied the common method (which comes with OpenCV samples), ie find contours in all color planes, approximate it and check for number of elements=4. It works to some extend that, few squares are detected, especially the dark ones.
Next step i did was prediction. ie this arrangement is fixed. So, if some are obtained, I can predict remaining ones. It also worked to some further extend. But accuracy was very bad.

But I feel prediction is not a good method here and it doesn't always provide accurate answers as given by the first step.

What I need :

1) Is there any other better methods to detect these squares more accurately? Or multiple methods?

One important point is that, time is not a problem here. Algorithm can be slow, it doesn't matter. But accuracy is the major criteria.

Sometimes, images can be much more blurry.

And one of the major problem I faced is that some squares have almost similar color as that of background (check column 3 first and second squares).

Looking for ideas, thanks in advance

UPDATE :

Below is the maximum accurate result I got :

Of course, the result image is resized a little bit.

UPDATE 2 :

I have given a much more better solution in my answer below: https://dsp.stackexchange.com/a/7526/818

A first attempt using Matlab:

im = imread('squares.jpg');
im2 = rgb2gray(im);

se = strel('disk', 15);

for i = 1:16;
    t = 60+i*5; % try out a range of bw thresholds to see what works best
    labelled = bwlabel(im2>t); % label regions in the BW image
    closed = imclose(labelled, se); % close small regions
    cleared = imclearborder(~closed,4); % clear regions touching the border
    subplot(4,4,i);
    imshow(cleared);
    title(['T = ' num2str(t)]);
end

Results in the following regions:

As you can see, selecting the threshold that results in the highest number of regions (T=120) would already give 7 correct locations, some merged locations, one false positive and two false negatives.

This was a fairly simple attempt but I think it shows that the approach works. Adding some stuff to break up elongated regions, or doing this for each color channel separately are just a couple of the things you could do to improve on this.

It would also help if you provided a few more test images.


gnuradio - flowgraph with cpp I build a flowgraph in gnuradio, generated python code, change a little and run it using python.
My question is, does building the flowgraph in cpp will provide better results (time)? (in other words does GIL and other stuff that make python slower also infect here?)(and if not i would like to read why)

Unless your flow graph contains signal processing blocks that you wrote in python (See the "embedded python block"), no, C++ will not make your flow graph be faster.
The reason is simple: Python is only used to set up the flow graph (read the generated python code!), i.e. to instantiate and connect the blocks together.
However, the GNU Radio scheduler and most blocks are running without interacting with the Python world at all. Hence, the GIL is needed at no point at all. Your Python calls top_block.run() and that's all that Python sees – the processing happens off-thread separately.


How does the OFDM system receiver decide the OFDM symbol length? I am trying to understand how the receiver decides the length of OFDM symbol in time domain before it goes into DFT process at the receiver.
The complete OFDM symbol (with Cyclic Prefix added) is propagated by the transceiver through the channel and it is received at the receiver with a time delay. This delay is caused by the propagation channel due to the objects or channel characteristics. And now lets say we have received the signal and first need to know what the symbol start time so that receiver can start the OFDM symbol and cut it at the end of the symbol. Then start for the consecutive symbol...

Could anyone please explain how to visualize this?

Thanks in advance!

First of all, the timing, phase and frequency offsets of the channel are removed. After synchronization, the cyclic prefix is removed, to avoid ISI between OFDM symbols. Now because the transmitter and receiver agree upon a common standard before the communication begins, the receiver knows how many samples to consider as a symbol.

Generally, the the number of samples in an OFDM symbol depend upon the bandwidth.


Distinguish vowels from consonants Problem of processing speech. Required to determine the phonemes and identify vowels and consonants. Anyone involved in this? Please advise what work on the subject is worth reading?

This is somehow a broad question. An interesting book is Digital Processing of speech signals by Rabiner & Schafer, in which you can find some good explanations to get started. It probably is out of date when looking into more advanced techniques, but I think it's good to understand the main characteristics of speech.

In short, every sound we produce is determined by the disposition of the vocal tract. In particular, if the vocal folds are tense and don't vibrate, the sound is said to be unvoiced. However, if they vibrate, the sound is voiced. You can sense the difference by placing your hand in your neck while you say "Ahhhh" or "Ssssshh". In Spanish, all vowels are voiced and all consonants are unvoiced (this is not totally true, please refer to the comments). I'm not sure about other languages though. You can think about a generator (glottis) which is connected to a filter (vocal tract) in a way that sounds take their characteristic shape. This is in fact the basis for speech production models (such as the simplest).

If you look closely to a speech signal, you can easily find differences between the different phonemes. Voiced phonemes are said to be pseudo-periodic. This means that, if you look at a fragment of the signal which is short enough (about tens of milliseconds long, e.g. 40ms), you can find the superposition of harmonic sine signals. This means that:

The signal is periodic with a fundamental frequency $f_0$. This is also known as the pitch.
Harmonics will appear with frequencies $2f_0, 3f_0, 4f_0$... which energy will depend on the exact phoneme.

If you plot the spectrum of this speech signal, you can find something like the following (source).

This resembles, for example, the harmonic structure of a brass instrument sound.

In contrast, unvoiced signals don't have such a clear harmonic structure. In fact, they have no fundamental frequency, and can be conceived as a noisy signal. You can sing by saying "Ahh", but you can't by saying "Ddd". Unvoiced sounds can be conceived as coloured noise. They usually have less energy than voiced.

Taking these differences into account, we can think about a way to detect them automatically.

  Voiced sounds are periodic while unvoiced aren't.

To detect periodicity, you can use methods such as autocorrelation (which has been widely treated in this site).

  Voiced sounds have more energy than unvoiced.

You can divide the signal in 20-40ms windows, compute the energy for each one and decide according to a threshold.

  Unvoiced sounds are noisy and tend to be prominent in high frequencies, while voiced aren't so.

Therefore, unvoiced will cross with zero more often than voiced. This is called Zero Crossing Rate.

While these simple methods are not too accurate by themselves, are a good place to start from. Furthermore, I am sure you can get good results by combining them. In any case, this post does not substitute a good book, as the one I mentioned before. It's just some kind of introduction to show you some possibilities.