texts: Wiener Khinchin theorem : struggle in the derivation

texts:
I am reading Introduction to quantum noise, measurement and amplification, and I need to understand the Wiener Khinchin theorem: how to derive it. I also need to understand some context around this theorem (why some object are defined the way they are).
The theorem is derived on the page 55 of this document. I will explain to you what I understand from the starting point of this theorem and where I am stuck.

I assume a stationary stochastic gaussian process. I call $R_{XX}(\tau)=\langle X(\tau) X(0) \rangle$ the autocorrelation function of a stochastic variable $X$. I assume the signals to be real.
We define the power spectral density of noise as:
$$ S_{XX}(\omega) \equiv \lim_{T \to +\infty} \langle |X_T(\omega)|^2 \rangle $$
Where $X_T(\omega)$ is the windowed Fourier transforme of $X$ on a period $T$. Basically we define the signal $X_T(t)$ that is equal to $X(t)$ for $t \in [-T/2,T/2]$ and that vanishes elsewhere. The Fourier transform of this signal can thus be defined as:
$$ X_T(\omega) = \frac{1}{\sqrt{T}} \int_{-T/2}^{+T/2} dt X_T(t) e^{i \omega t} $$
To be consistant we must define the inverse Fourier transform as:
$$ X_T(t) = \frac{\sqrt{T}}{2 \pi} \int_{-\infty}^{+\infty} d \omega X_T(\omega) e^{-i \omega t} $$
First questions:

Why do we take such "uncommon" definition for the Fourier transform
(usually it is $\frac{1}{\sqrt{2 \pi}}$ for both sides
(direct+reciprocal), or $\frac{1}{2 \pi}$ in front of one of the
sides. Is there a theoretical motivation for this in this context ?
Why do we need to window the function ? Is it to avoid technical difficulties with divergent integrals ?

Now, I go to the derivation and express where I am stuck. The goal is to show that $R_{XX}(\omega)$, the Fourier transform of the autocorrelation function corresponds to the noise power spectral density $S_{XX}(\omega)$. I define $S^T_{XX}(\omega)$ the quantity before the limit of $T \rightarrow +\infty$. I know how we can derive computing the inverse Fourier transform of $S_{XX}(\omega)$ directly as suggested in a link in the comment. However I am stuck on deriving it from the other side (I have $S_{XX}(\omega)$ on the left hand side and I make appear $T.F(R_{XX})$ on the right handside). It might look as a boring question but in the document I follow they always follow this way of proving equalities, and I really don't get the manipulation they do. To understand how they do I need to understand this "other way" of deriving the theorem.
$$ S^T_{XX}(\omega)=\frac{1}{T} \int_{-T/2}^{+T/2} dt \int_{-T/2}^{+T/2} dt' \langle X_T(t) X_T(t') \rangle e^{i \omega (t-t')}$$
Now, I use the fact my process is stationnary. And I recognize the autocorrelation function:
I start from the definition (before taking the limit):
$$ S^T_{XX}(\omega)=\frac{1}{T} \int_{-T/2}^{+T/2} dt \int_{-T/2}^{+T/2} dt' R_{X_T X_T}(t-t') e^{i \omega (t-t')}$$
I do the change of variables: $(t,t') \rightarrow (\tau=t, \tau'=t-t')$. It gives me:
$$ S^T_{XX}(\omega)=\frac{1}{T} \int_{-T/2}^{+T/2} d \tau \int_{t-T/2}^{t+T/2} d \tau' R_{X_T X_T}(\tau) e^{i \omega \tau}$$
At this point, I have absolutely no idea how we could proceed.
We can also try another method (which I think is the way they do in the document but I really don't understand their justification). We can do the change of variable $(t,t') \rightarrow (v=t+t',u=t-t')$. The Jacobian gives $1/2$ and the boundaries are the one I write below:
$$S^T_{XX}(\omega)=\frac{1}{T} \int_{-T}^{T} dv \int_{-B(v)}^{+B(v)} \frac{du}{2} R_{X_T X_T}(u) e^{i \omega u}$$
$B(v)=v+T$ for $-T \leq v \leq 0$ and $B(v)=t-v$ for $0 \leq v \leq T$. It gives:
$$S^T_{XX}(\omega)=\frac{1}{T} \int_{0}^{T} dv \int_{-(T-v)}^{+(T-v)} \frac{du}{2} R_{X_T X_T}(u) e^{i \omega u}+\frac{1}{T} \int_{-T}^{0} dv \int_{-(T+v)}^{+(T+v)} \frac{du}{2} R_{X_T X_T}(u) e^{i \omega u}$$
And same problem here: I really don't see how we could get to the result with this way of deriving (which I think is roughly the path taken is the document I linked)
Thus, second question: How can prove the theorem with this method ?


texts:
I'm trying to understand DFT in the context of an integer ring (i.e. the number theoretic transform) and I was wondering how to reduce the number of multiplications by the twiddle factors.
Consider 3-point DFT:
X_0 = x_0 + x_1       + x_2
X_1 = x_0 + x_1 * w   + x_2 * w^2
X_2 = x_0 + x_1 * w^2 + x_2 * w^4

Which is the same as:
X_0 = x_0 + x_1       + x_2
X_1 = x_0 + x_1 * w   + x_2 * w^2
X_2 = x_0 + x_1 * w^2 + x_2 * w

and it needs 4 (modular) multiplications.
However, apparently it can be also written as:
X_0 = x_0 + x_1           + x_2
X_1 = x_0 + x_1 * w       - x_2 - x_2 * w
X_2 = x_0 - x_1 - x_1 * w + x_2 * w

which uses only 2 multiplications.
I just cannot figure out how to derive - x_2 - x_2 * w out of x_2 * w^2, could someone please fill in the details for me? Are there other ways to reduce the number of multiplications? Is it possible to use just one (or even zero) multiplication? Does similar trick work for larger transforms?


texts:
I‘m a high school student and I haven’t studied physics or anything.

Why does the DFT depend on an integer, say $k$ or $n$ (it’s usually expressed like $F(n)=...$ or $F(k)$ or $F_k$, etc.) if it is supposed to deliver a frequency information of a sampled signal?

Can the frequency content of the signal be expressed as a multiple of the integer?


texts:
I've implemented Goertzel algorithm, now very simply in Python. But I cannot obtain the correct answer it's supposed to produce for a single DTF frequency. here is my code in Python:

import numpy as np

# omega frequency: radians/sample
omega = 0.23

# length of the data to generate
length = 200

# updated parameters
dft = 0.0
s1 = 0.0
s2 = 0.0

# initial phase
radians = 0.0
for i in range(length):

    # generate signal value
    value = np.cos(radians)

    # goertzel algorithm update
    s0 = value + 2.0 * np.cos(omega) * s1  - s2
    s2 = s1
    s1 = s0

    # dft sum update
    dft += value * np.exp(-1j * omega * i)

    # next value
    radians += omega

# final goertzel algorithm update
s0 = 2.0 * s1 * np.cos(omega) - s2
goe = s0 - np.exp(-1j * omega) * s1

# print results
print "DFT: {:.4f}".format(dft / length)
print "GOE: {:.4f}".format(goe / length)

Output:

DFT: 0.4979-0.0097j
GOE: -0.2065+0.4531j

The correct answer (with infinite data) is 0.5 + 0.0j.

DFT is therefore providing the correct result.

What's wrong with my Goertzel algorithm implementation?

https://en.wikipedia.org/wiki/Goertzel_algorithm


texts: Can sample aliasing of strictly white noise lead to random walk (drift) in the sampled signal?

texts:
Aliasing of tones higher than the Nyquist frequency, and where the tone may have some degree of nonlinear distortion can appear as low frequency 'drift' or even near static bias in the sampled signal.

Also it's known that integration of a signal with a large component of white noise results in random walk - a 'drift' like behavior in the integrated signal.

So are there any possible situations combined with aliased sampling of signals with a significant white noise, or pink noise that can appear as low frequency drift in the sampled signal?


texts: Mathematical relationship between highpass and lowpass filtering

texts:
Let $g, h_{HP}, h_{LP}: \mathbb{R} \rightarrow \mathbb{R}$ and $G, H_{HP}, H_{LP}$ denote their continuous Fourier transforms under the Fourier operator $\mathcal{F}$. Let $*$ denote the continuous convolution of two functions and $\cdot$ the pointwise multiplication.

Define the shorthand $1: \mathbb{R} \rightarrow \mathbb{R}, x \mapsto 1$ and let $H_{HP} = 1 - H_{LP}$ in the frequency domain, then I would like to show
$$
g * h_{HP} = g - (g * h_{LP}).
$$

In words:

  Filtering an image with a highpass is the same as subtracting the lowpass-filtered image, where the lowpass is the inverted highpass.

There is a proof (from some closed-source course notes) which goes as follows:

\begin{align*}
& g * h_{HP}\\
& = \mathcal{F}^{-1}(G \cdot H_{HP})\\
& = \mathcal{F}^{-1}(G \cdot (1 - H_{LP}))\\
& = g * \mathcal{F}^{-1}(1 - H_{LP}) \quad \text{(The next step is unclear to me)}\\
& = g * (1 - h_{LP})\\
& = g - (g * h_{LP})
\end{align*}

In fact, $\mathcal{F}^{-1}(1 - H_{LP}) = 1 - h_{LP}$ is equivalent to $1 = \mathcal{F}(1)$ by applying $\mathcal{F}$ on both sides, which is a contradiction with the first identity mentioned below.

Questions

  Is it true that the proof above lacks details?
  Is my proof below correct/more rigorous?

It is known that

$\mathcal{F}(1)=\delta$, from Constant function on MathWorld
$\mathcal{F}^{-1}(1)=\delta$, from Math Stack Exchange: How to prove that inverse Fourier transform of “1” is delta funstion?

I found another approach to possibly prove the proposition above:

First note: If $u, v: \mathbb{R} \rightarrow \mathbb{R}$ are continuous functions and $\int_a^b u(x) dx = \int_a^b v(x) dx$ for every $a, b \in \mathbb{R}$, then $u(x) = v(x)$ for every $x$, i.e. $u=v$.

\begin{align}
&\int_{a}^{b} (g * h_{HP})(x) dx\\
& = \int_{a}^{b} (g * \mathcal{F}^{-1}(1 - H_{LP}))(x) dx \quad \text{(as stated above)}\\
&= \int_{a}^{b} (g * (\delta - h_{LP}))(x) dx \quad \text{(linearity + identity above)}\\
& = \int_{a}^{b} \int_{-\infty}^{\infty} g(\tau) ((\delta - h_{LP})(x - \tau)) d\tau dx\\
& = \int_{a}^{b} \int_{-\infty}^{\infty} g(\tau) (\delta(x - \tau) - h_{LP}(x - \tau)) d\tau dx\\
& = \int_{a}^{b} \int_{-\infty}^{\infty} g(\tau)\delta(x - \tau) - g(\tau)h_{LP}(x - \tau) d\tau dx\\
& = \int_{a}^{b} g(x) - \int_{-\infty}^{\infty} g(\tau)h_{LP}(x - \tau) d\tau dx\\
& = \int_{a}^{b} g(x) - (g * h_{LP})(x) dx\\
& = \int_{a}^{b} (g - (g * h_{LP}))(x) dx
\end{align}

Especially, I am unsure whether the usage of the $\delta$ distribution is rigorous here.


texts: Filtering an image with two gabor filters

texts:
I have an image of size 160x160 of pixels. I am trying to convolve this image with two Gabor filters of different scales (7 and 9 respectively) but with the same orientation. I red in an article the following phrase :

The maps (the results of the two convolution operations) have the same dimensionality but they are the outputs of different filters.

After convolving my image (conv2) with these two gabor filters (scales 7 and 9), i got the following results :

The 2 maps (one image after convolving with gabor of scale 7 and the other after convolving with gabor of scale 9) have different dimensionality!

So according to the phrase 1) , why my maps have different sizes ? it is normal ? or false ? i need your opinions.

In other words, i filtered my image with a gabor filter of scale 7, then i filtered my image with another gabor with scale 9 . So i got two maps. Is these maps must have the same size or not ?


texts: Problem with Covariance matrix using diagonal loading involved in calculation of eigenvalues

texts:
I have the following problem : I'm calculating the sample covariance matrix in the frequency domain ( $y_{k}$ is the FFT of a time domain $k_{th}$ symbol  vector signal , basically a simulated received signal) as follows:

$$\mathbf{R}=\frac{1}{N_{f}}\sum_{k=0}^{N_{f}-1}y_{k}y_{k}^{H}$$

Well , the next step on my algorithm is to solve a optimization process,  as an essential part of it I need to compute the eigenvalues by doing SVD ,method of powers, etc... in MATLAB of the following equation:

$$(\mathbf{R}^{-1}\mathbf{A})$$

To not get into too many details because I believe my issue comes from a numerical problem and many insights of the actual algortihm / context aren't needed .Let's assume, $A$ is simply a predefined matrix that I compute.

The REAL ISSUE appears now because $\mathbf R$ is ill-conditioned as MATLAB tells me so. So the inverse procedure seems to be failing and the eigenvalues I'm obtaining are incredibly small due to this issue (In fact I only get 1 eigenvalue different than zero). The dimensions of $\mathbf R$ are typically large ( since they are compressed depends on the actual compression ratio I'm using but let's say $32\times 32$ for example).

One approach to solve this problem I found is to use diagonal loading as:

$$(\left(\mathbf R+\sigma\mathbf I\right)^{-1}\mathbf A)\quad\text{with}\quad\sigma > 0$$

This seems to be solving the problem, the eigenvalues are now scaled due to this background "noise". My question is how can I obtain the truly well scaled eigenvalues because , later on,  in my algorithm these eigenvalues will serve as weights since I'm considering them as an actual power estimate.

Note : I've been playing with the cond() function in MATLAB for $\sigma= 0.05$ ,cond(R+sigma*I)= 2 , which is not bad I believe.

Feel free to ask more questions about the problem. But I think my question relates to a purely numerical issue involving eigenvalues and covariances matrices.


texts:
My colleague and I started exploring EEG processing and, especially, functional connectivity. We decided to try the Phase Slope Index (PSI) and conducted a simple experiment as follows: it consists of three steps, research participant should raise his left hand (1), then his right one (2), then both (3). The description is simplified a bit but shows what the experiment was. We expected to see the information flow between the pre-motor cortex and primary motor cortex but we seemed to fail. The data was pre-processed (filtered at [0.25; 40] Hz, re-referenced to the average electrode, segmented into three pieces according to three stages of experiment; also blink artefacts were removed, so were motion ones). Then we computed the all-to-all PSI matrix and what we have had is depicted on the pictures (if needed and more appropriate, I could add the graph representation).   All in all, we haven't seemed to reach what we aimed and now we're trying to figure out whether the fail refers to the interpretation, or experiment design, or both, or something else. My questions are:

Can one guess what stage of experiment every matrix refers to? If so, how? I mean, do they at least look like the result of the experiment described? If they do, how can one interpret them?
What could be an experiment that would show the noticeable PSI for two areas given beforehand? I mean, we had two areas given, pre-motor cortex and primary motor cortex, and we supposed that we would find a synchronization between them in such a motion experiment; something's gone wrong, though. If the problem is in the design of experiment, what design would be better?


texts:
This might be a silly question. I was reading about the Laplacian of Gaussian (LoG) operator and got confused about the alternative equivelant ways we can make use of it.
Let's assume we have a 2D image $I(x,y)$
Since linear convolution commutes with differential operators the following equality is valid:
$$\nabla^2 (G_{\sigma}\ast I) = (\nabla^2 G_{\sigma})\ast I $$
We also know that convolution is commutative, so we could also write (could we??) :
$$ \nabla^2 (G_{\sigma}\ast I) = \nabla^2 (I \ast G_{\sigma})  \mathop{=}\limits^?G_{\sigma}\ast(\nabla^2 I)$$
which means that we first have to  apply the Laplacian operator on image $I(x,y)$ and then convolve the result with the Gaussian filter.
However, is this approach correct? It seems to me wrong and haven't found it anywhere. I tend to believe it's not valid since it will amplify the noise of initial image $I(x,y)$ (because we first apply the Laplacian operator on the noisy image) and then the smoothing will not be such efficient.


texts: Differences between DSP programming vs MPU programming

texts:
I'm moderately skilled in programming classic microcontrollers. Embedded C, built-in peripherals, interrupts, the daily bread of an Embedded programmer.

I'm interested in expanding my skillset onto programming DSP. Of course the basic steps would be grabbing some development kit, some manuals and tutorials and get to writing.

But I've seen microcontroller code written by a seasoned Java developer starting with embedded. I've seen an Embedded engineer trying to write web apps. I myself moved from desktop applications development to writing backends running on server farms with thousands requests per second.

The results are not pretty - the paradigms, habits, styles, hang-ups result in using design patterns simply not applicable, optimizing elements that given platform should not have optimized, and so on. It took some tutoring from someone experienced to learn the necessities and unlearn the wrong parts.

So - as a beginning programmer of DSP with background in standard MPUs (ARM7, 8052 etc) what must I learn, and what must I unlearn to write good DSP code? What are the basic differences between the styles and prerequisites of programming the two?


texts: Me and a coworker are currently in the process of making an audio dsp library. We are wondering what some good FIR and IIR filters to use were. I know ALOT of them but what do you guys think in best/necessary for audio? Also if you have any other ideas not related to FIR and IIR filters please give your input.


texts: make a threshold frequency based on percentage of frequencies less than that?

texts:
Using MATLAB signal processing toolbox, I have created signals like dirichlet, chirp and etc and and obtained their fourier and cosine fourier transform. Now I want to find the frequency where for example A% frequencies are less than that and use it as a thresholding frequency to compact the signal and eliminate unnecessary details?
Could you please tell me a little more about thresholding frequencies and how to choose them in signal processing?


texts: Creating realization of 2D Gaussian field in Fourier space

texts:
I want to generate a 2D Gaussian field with dimensions $L\times L$ with $N^2$ cells each of size $l = L/N$. I'm doing this by producing a realization of this field in Fourier space by producing Gaussian random numbers with variance $\sigma^2(k) \propto P(k)$ where $P(k)$ is some input power spectrum.

To do this I need to have a grid in Fourier space with each cell defining a $k$-value. But I'm not sure how each component $k$-value combines in each cell, i.e. each cell will presumably have a $k_x, k_y, k_z$ value depending on the cells position in Fourier space.

Do I combine them by doing $\textbf{k} = \sqrt{k_x^2 + k_y^2 + k_z^2}$? But then surely $\textbf{k}_{max}$ will then be too large in the corners i.e. larger than $2\pi/l$

Do I define different values of $\sigma$ i.e. $\sigma_i(k_i) \propto P(k_i)$ for different $x,y$ and $z$ then combine all these different $\sigma$'s when generating Gaussian random numbers.

I'm also confused by the range of values that my $k$ should have since examples such as the attached image suggest that I should include a value of $k=0$ (representing the DC frequency) but I can't then generate a Gaussian random number for this cell since my power spectrum will be undefined at $k=0$.

Layout of discrete Fourier space


texts: Why do we always use square kernels for filters?

texts:
In image processing, why do we always use odd square kernels(3x3, 5x5, 7x7, etc.) for filters?

Why can't we use kernels of size 3x1, 5x3, etc. that are rectangular kernels? Also, why do we not prefer even kernels (2x2, 4x4, 6x6, etc.)?


texts: Moving average on the frequency domain - (Welch's Method?)

texts:
Is there a name for a filter that is using a moving average on the frequency domain?

I am searching for such a filter to filter out long "pads" (synths) from music.

Do any of you know about such software?


texts:
I want to implement this function $[n,m] = a^{n+m}u[n]$ using 41×41 array, where $u[n] =1$ since all the value will be positive, and $n,m$ from 0 to 40. But it gives a wrong value.

N=40
M=40
a=0.9
x =zeros(41,41)
for col = 1:41
  for row = 1:41
    for k=0:N
      for i=0:M
        x(row,col)=a.^(k+i)*1;
      end
    end
  end
end


texts: Preprocess sound for voice recognition in Pocket Sphinx on Raspberry Pi 3

texts:
I want to improve my voice recognition system on a Raspberry Pi 3 running Pocket Spinx. It's pretty accurate (limited recognition vocabulary) but is terrible when other noise is around.

Is there any way to pre-process sounds before feeding them into pocket sphinx to filter only the human voice ? I was thinking there might be some way to reduce the data being passed in, eg human voice being between certain frequency etc, white noise removal etc. I haven't found anything simple so far. I'm thinking of filtering and processing the voice recordings on an STM32 that has an Arm Cortex M4 and DSP etc.


texts:
Sorry if this sounds elementary but I am struggling to grasp the physical idea behind ARMA (auto-regressive, moving average) process. The "AR" part is intuitive and so is "MA", but put together?

If I model my time series using "AR", I can predict the next sample using linear terms of the previous samples and the difference between the predicted and the real one should then have a normal distribution provided the noise was guassian white (the unpredictable).
But what does it mean to say that my time series can be ARMA modelled?

thanks


texts: Convolution by differentiation property of convolution

texts:
This is the question

when i am trying to do it by the deravative property of convolution
https://en.wikipedia.org/wiki/Convolution#Differentiation like this

Actual answer should come like this

  $$y(t)=
\begin{cases}0, & \text{for }t < 0 \\
7-3t, & \text{for } 0<t<1 \\
5-t, & \text{for } 1<t<5 \\
7, & \text{for } t>5
\end{cases}$$

i am getting answer -6 at t=4 bu the actual answer should come 1

  P.S:-
  Please dont give me alternative method to solve this i know it already instead tell me whats wrong i am doing in this method
  also i know i am violating the policy by putting the image but i will edit it lator in text form .


texts: Determining Loss of Information by Taking Average (Mean of Signal)

texts:
So basically information is defined by expected value of Shannon's information i.e. Entropy. I am curious how much information is lost if we simply take the average of the sample given to us. I am also curious if we also take the second moment (variance) of the same signal, how much information do we gain by using the variance along with the mean? In other words:

$$H_{\textrm{Mean}}(\Sigma_0^n x_iP_i)=\ ?$$
$$H_{\textrm{Variance}}(\Sigma_0^n (x_i-\mu )^2)=\ ?$$

Motivation: I just want to know how much information is dropped if we simply drop variace term and use just mean instead. I believe this concept cant be expanded to also consider information in third and fourth moments.


texts:
The $\frac{d_1}{d_2}<0.6$ gives a large number of false matches and it assumes the nearest neighbour is a correct match.

Anybody knows any good method to reduce number of false matches or filter them ?

Can we use k-NN matches instead of $\frac{d_1}{d_2}<0.6$ ratio for matching the SIFT features and then choose the best of k points by:

Changing the Euclidean distance by considering the spatial neighbourhood, or
Spatial neighbourhood matching.

Anybody ever tried something like that? Or any paper related?


texts:
STFT can be successfully used on sound data (with a .wav soundfile for example) in order to do some frequency-domain modifications (example : noise removal).
With N=441000 (i.e. 10 seconds at sampling rate fs=44100), windowsize=4096, overlap=4, STFT produces approximatively a 430x4096 array (first coordinate : time frame, second coordinate : frequency bin). Modifications can be done on this array, and reconstruction can be done with overlap-add (*).

How is it possible to do a similar thing with wavelets ? (DWT), i.e. get a similar  array  of shape a x b, with a time frames, and b frequency bins,  do some modification on this array, and at the end, recover a signal ? How ? What is the wavelet equivalent to overlap-add ?
What would be the Python functions involved here (I haven't found an easy example of audio modification with pyWavelets...) ?

(*) : Here is the STFT framework that can be used :

signal = stft.Stft(x, 4096, 4)    # x is the input
modified_signal = np.zeros(signal.shape, dtype=np.complex)

for i in xrange(signal.shape[0]):    # Process each STFT frame
    modified_signal[i, :] =  signal[i, :] * .....  # here do something in order to
                                                   # modify the signal in frequency domain !
y = stft.OverlapAdd(modified_signal, 4)   # y is the output

The goal is to find a similar framework with wavelets.


texts:
I have a question regarding the inverse Fourier transform and its relevance to non-stationary signals. And by non-stationary signal, I'm talking about a signals whose frequency content varies with time. I'll provide the example that made me originally ask this question.

Say for example I have a sinusoidal signal with a single frequency that goes on for ever. I could take the Fourier transform of this signal to obtain its frequency content. And because this signal is a single frequency that goes on forever, this signal is considered stationary.

However, now let's say I have a signal that is zero except for in some window in which the signal is sinusoidal. This nonzero window exists from t = 0 to 1 second. If I were to take the Fourier Transform of this signal over a finite duration (such that it acts on both the nonzero and zero part of the signal), I would get a frequency content that shows there exists a peak frequency content due to the nonzero window.

Now let's say I were to do the same thing, except this time the nonzero window exists from t = 6 to 7 seconds. If I were to take the Fourier transform, over the same finite duration, I believe I would yield the exact same frequency content (both phase and magnitude) from the earlier experiment.

Assuming what I said was correct, it should thus be IMPOSSIBLE to recover the exact signal from the Fourier transform since two different time domain signals have mapped to the same frequency domain signal. Am I correct with this statement? Furthermore, I believe the inverse Fourier transform will end up producing a stationary signal. As in, the frequency that was once contained just within the nonzero window will now be shown to be infinite from the inverse Fourier transform.

TLDR:

Does the inverse Fourier transform only produce stationary time signals?
Thus if I have non-stationary signal, take the Fourier transform and then take the inverse Fourier transform, I will now have a signal that is stationary?


texts: Normalized cross-correlation in frequency domain

texts:
I never worked with signal processing and never really used Fourier transforms before, still I am working on a project consisting on taking the output of an accelerometer to detect some movement features between touching surfaces. Using a frame of the output, I compute the FFT. After analyzing the data for multiple experiments, it appeared to make sense to use gaussians to isolate bands of frequencies and use that to classify the features.

That does not work well, however, if I change the sensor I am using, the sampling rate or frame size. Thus, I want to calculate the normalized cross-correlation. Problem is: I could not find an algorithm to calculate the NCC in the frequency domain and I did not figure out how to transform the translated gaussians to the time domain.

For instance, I have a frequency spectrum from -F/2 to F/2 and L samples. I converted that to L/2+1 frequency bins from 0 to F/2 and defined a band using a gaussian with standard deviation s and centered on c.
To get the cross-correlation between the spectrum and the gaussian, I simply do the sum of the element-wise multiplication of them. Now, if I want that normalized, what can I do? Is there a way to normalize it right from the frequency domain?


texts: Formulating a function on Matlab for the Shannon interpolation formula

texts:
I am trying to formulate an algorithm for applying to Shannon interpolation formula to the discrete signal

$$x[n]=\frac{c^2}{4}\int\limits_{0}^{n T}y\left(\tfrac{c}{2}s\right)\,ds,$$

where $c$ is constant. Now, if $\frac{1}{T}=f_{s}>2B$, where $B$ is the band-limit of $X(f)$, then we can reconstruct the continuous time signal as:

$$x_{\text{const}}(t)=\sum\limits_{n=-\infty}^{\infty}\frac{c^2}{4}{\int\limits_{0}^{n T}y\left(\tfrac{c}{2}s\right)\,ds}\cdot\text{sinc}\left(\frac{t-n T}{T}\right)$$

function x=interpolate(y,N,T,c)
     x=0;
     for n=1:N
         x=x+((c^2)/2)*integral(y,0,n.*T).*sinc((t-n.*T)./T);
     end

However, I am having trouble with incorporating the t variable here. Setting it as global and running interpolate(y,20,20,1) also produces nothing to note.

Note that the continuous time signal in question is given by

$$x_{\text{true}}(t)=\frac{c^2}{4}\int_{0}^{t}y\left(\tfrac{c}{2}s\right)\,ds.$$


texts: Normalized cross correlation via FHT - how can I get correlation score?

texts:
I'm using the 2D Fast Hartley Transform to do fast correlation of two images in the frequency domain, which is the equivalent of NCC (normalized cross correlation) in the spatial domain.

However, with NCC, I can get a confidence metric that gives me an idea of how strong the correlation is at a certain offset. In the frequency domain version, I end up with a peak-finding problem in the inverse FHT after doing the correlation, so my question is:

Can I use the value of the peak that I find in the correlation image to derive the same (or similar) confidence metric that I can get from NCC? If so, how do I calculate it?


texts:
I'm trying to perform a FFT of a 3D regular grid and then compute the bin average (in spherical shell bins) of the Fourier transformed grid.

The problem is that the resulted vector is very noisy as I'm averaging over just a few frequencies in each bin.

To correct this and get a smoother outcome, I tried zero-padding the initial grid before transforming since I've read that this equivalent to Sinc interpolation.

Now the problem is that I'm running out of memory.

Is there a way to perform the interpolation (Sinc or not) without zero padding? Let me also point out that this is a 3D FFT.

I'm doing the whole thing in Python.

Update

Let me try to simplify the question a bit:

You can ignore the spherical shell bin-averaging part.

I know that using zero-padding can increase the number of frequencies that you can get. However, this comes at the cost of increased memory usage. Due to the fact that I'm using 3D FFT (not just 1D) memory is a real issue.

So my question is: Can I have the benefits of zero-padding (many frequencies) without its computational burden by using some other method? (e.g. interpolation) And if so, then how?


texts: Adding noise in frequency domain while solving an SPDE

texts:
I’m implementing a numerical solver to a certain kind of a 2D stochastic PDE using Fourier Transforms.
So, this solver is a Galerkin type method, which means that choosing the right basis functions I can apply Fourier Transforms to compute the solution.
Long story short, applying the explicit Euler-Maruyama scheme to my SPDE and discretising the domain with $N^2$ nodes, the discrete version of the Fourier coefficients of $V(t_{i+1})$ is given by:
\begin{equation}\tag{1}
\hat{v}_{mn}\left(t_{i+1}\right) = \hat{v}_{mn}\left(t_{i}\right) + \frac{\Delta t}{\alpha}\Big(\hat{\imath}_{mn}\left(t_{i}\right) - \hat{v}_{mn}\left(t_{i}\right) + \hat{a}_{mn}(t_i)\Big) + \epsilon \sqrt{\Delta t} \lambda_{mn}W\left(t_i\right),
\end{equation}
where the hat terms corresponds to the Fourier Coefficients of the equation's functions, $\epsilon$ is the level of additive noise, $W$ is a complex random matrix with distribution $N(0,1)$ and $\lambda$ is a spatial correlation matrix.
Notice that the noise is being added in the frequency domain!
However, when comparing results with other numerical methods the effect of the noise is not the same, actually I have to scale up $\epsilon$ by a factor.
I'm using the fftw library in julia to compute the Fast Fourier Transforms and it is known that this library do gives us the Fourier coefficients multiplied by $N^2$, i.e. $fft(V) = N^2 \hat{v}_{mn}$. So, the stochastic term of (1) has not the same magnitude as the rest of the coefficients. Am I right?
To overcome this issue, I simply have to multiply $\epsilon$ by $N^2$, is this easy? Or adding noise in the frequency domain is more complicated than this?
Sorry for my question, but I’m a bit newbie in the Fourier world.
Actually, an experimental evidence that the number of nodes influences the noise effect in the solution, is if I increase $N$ I'll also have to increase $\epsilon$ to continue to have similar behaviours than with less nodes.
Thank you!
Tiago


texts: Temporally synchronize two different or more real time digital audio signals

texts: Are there any DSP techniques to temporally synchronize two different or more real-time digital audio signals without timestamp information?


texts: Why the radar (moving) target return is not centered at zero frequency

texts:
In Radar signal processing, one has to separate clutter(unwanted signals) from moving target signals.
The picture given below is an example of Radar Signal Power Spectum(taken from the book Radar Systems Analysis  and Design Using MATLAB by B.R. Mahafza).

Why target return is coming in between 0(DC) and $f_{r}$ where $f_{r}$ is the pulse repetition frequency.


texts: How to properly match bit growth of a FIR filter with a reference value

texts:
I am trying to implement a chain of CIC/FIR filters on an ZYNQ FPGA.
Using the Xilinx FIR compiler works fine so far but I am unable to properly get all the math.

At the moment I have 2 chains of filters:

$5\times$ downsampling, passband $12.5\text{MHz}$:
→ FIR filter 2, 5x downsampling, input 17.7.0, output 17.7
$625\times$ downsampling, passband $100\text{kHz}$:

→ CIC filter, $25\times$ downsampling, input 14 bit, output 16 bit

→ Compensation FIR filter, no downsampling, input 17.7, output 17.7

→ FIR filter 1, $5\times$ downsampling, input 17.7.0, output 17.7

→ FIR filter 2, $5\times$ downsampling, input 17.7.0, output 17.7

I have a 14 bit, $20\text{kHz}$ signal that I feed into both chains:

If I feed it into the first chain I sign extend the 14bit value to 17bits and add 7 zero-bits at the end (.7).
I get a nice output that wins me one bit.
If I feed the signal into the second chain I directly feed the 14bit signal into the CIC filter and left-shift the output signal to only keep the 16 MSBs. I then sign extend that 16bit value to 17bit (yes I could have also just kept the first 17bits of the CIC return value, if that would be better, please tell me) and add 7 zero-bits at the end (.7) as well.
The filters win me 1-2 bit until the output of the "FIR filter 1" but at the output of the "FIR filter 2" I have an output where I lost 8bit (upper 8 non-sign bits are 0). This is the exact same FIR filter that gives me +1 bit in the first ($5\times$ downsampling) chain ... This broken chain also seems to have a lot of jitter after one period ... (I can exclude overflows as the source of error as the Xilinx FIR compiler guarantees no overflows).

How can this be, what am I missing?

In the same question I would like to ask:

Let's assume all the chains finally work properly. I then will have different bit growth. Initially I knew that 14bit correspond to $2.2\text{Vpp}$ which means that 12.86bits correspond to $1\text{V}$. But now I suddenly have a 17.7bit output. How do I know where my $1\text{V}$ reference level is now? Because it has shifted some.
I am quite sure that accounting just for the DC gain yields wrong results.

My waveforms:

Output $5\times$ chain:

Chain $625\times$:

I am very grateful for any hints and will gladly add more information if some more is required.


texts: Relationship between $E_b/N_0$ and $\rm SNR$ with channel coding

texts:
I would like to know what would be the relationship between $E_b/N_0$ and $\rm SNR$ when there is a coding rate of $R_c=K/N$. Normally when there is no channel coding according to this link.

$$\textrm{SNR(dB)} = \frac{E_b}{N_0}(\textrm{dB}) + 10\log10\left(\frac{R_b}{B}\right)$$

However what would happen if I have a linear block code or other coding schemes with code rate $R_c =1/2$ or other values, how to calculate the relationship between $\rm SNR$ and $E_b/N_0$ in this case ?


texts:
For quality control, I want to measure the frequency response of hardware I built. What I can easily generate on the hardware is a step response. Thus, I want to derive the magnitude and phase from the step response. This works well in general, however I am looking for ways to improve and understand basic relationships.
Remark: Quality control is probably done best by comparing the measured step response directly with a "golden" step response. However, I want to provide also numbers and visuals of magnitude and phase as the original system specification is given in these terms.
Find below the python code of what I have, comparing my approach with a direct Bode calculation from the transfer function.
Questions:

Is there a formula or rule of thumb of how the duration/length in time of the recorded step response and the sample rate used affect the quality/uncertainty of magnitude and phase? One obvious thing is, that Bode data for low frequencies with longer period than the recorded duration will have lower quality.
I always get a phase jump at the nyquist frequency, regardless of the sample rate and signal duration used. I am struggling to identify the root cause of this. Does it have an explanation from DSP theory or is it my implementation?
What other things should I pay attention to to get high quality Bode data?

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import signal # to generate example data

def stepResp2Bode(stepResponse, fS, minFreqRes=None):
    if minFreqRes is not None:
        if fS/stepResponse.size > minFreqRes:
            print('Minimum frequency resolution not achieved!')
            return # not enough samples for required frequency resolution
    impResp = np.diff(stepResponse, axis=0, prepend=0) # get impulse response by differentiating step response
    fftResult = np.fft.rfft(impResp, axis=0) # do fft
    bodeF = np.fft.rfftfreq(impResp.size, d=1/fS) # create frequency vector for fft-Results
    bodeAmplDB = 20 * np.log10(np.abs(fftResult)) # extract amplitude in dB
    bodePhase = np.angle(fftResult) # extract phase in radians
    results = pd.DataFrame({'frequency':bodeF,
                        'amplitude in dB':bodeAmplDB,
                        'phase in radians':bodePhase,
                        })
    return results

if __name__ == "__main__":
    # generate example data - in the application this would be real world measured data
    fS = 50 # sampling frequency
    wc = (2*np.pi*1) # corner frequency of example low pass
    t = np.arange(0, 5, 1/fS) # vector of time values
    sys = signal.lti([1.0], [1/wc, 1.0]) # generate lowpass filter
    stepResponse = signal.step(sys, T=t)[1] # calculate step response of low pass filter at given time values
    # Calculation
    bodeStepResults = stepResp2Bode(stepResponse, fS, minFreqRes=1) # from step response to magnitude and phase
    w, mag, phase = signal.bode(sys) # bode plot directly from system description
    bodeScipyResults = pd.DataFrame({'frequency':(w/(2*np.pi)),
                    'amplitude in dB':mag,
                    'phase in radians':(phase/180*np.pi),
                    })
    # plot results
    fig = plt.figure(1)
    fig.clear()
    ax1 = fig.add_subplot(211)
    ax1.semilogx(bodeStepResults['frequency'], bodeStepResults['amplitude in dB'], '-', label='from step response')
    ax1.semilogx(bodeScipyResults['frequency'], bodeScipyResults['amplitude in dB'], '-', label='from signal.bode')
    ax1.set_xlabel('Frequency  in Hz')
    ax1.set_ylabel('Amplitude  in dB')
    ax1.grid(which='both', axis='both')
    ax1.legend(loc='upper center')
    ax2 = fig.add_subplot(212)
    ax2.semilogx(bodeStepResults['frequency'], bodeStepResults['phase in radians'], '-', label='from step response')
    ax2.semilogx(bodeScipyResults['frequency'], bodeScipyResults['phase in radians'], '-', label='from signal.bode')
    ax2.set_xlabel('Frequency  in Hz')
    ax2.set_ylabel('Phase in radians')
    ax2.grid(which='both', axis='both')
    ax2.legend(loc='upper center')
```


texts: I'm an Electrical Engineer working as an independent consultant on applied research type DSP projects at small startups. Can anyone please recommend a good example contract to use with my clients? Are there any good templates out there which balance the rights and responsibilities of both parties?


texts: Is it possible to design a digital filter that rejects 60Hz noise while keeping step response in the time domain below one mili-second

texts: Have been trying to design IIR and FIR 60Hz notch filters but all have step response way above 10ms.


texts: how to lower frequency of a sound file in matlab?

texts:
I am new to programming and sound processing,
I have been using matlab to process sound files, a few of which are ultrasonic.
my sampling rate is 250kHz.
I Would like to be able to play the sounds at an audible frequency, there for i need to lower the pitch of the sounds without changing the speed.
The frequency in the sound files are around 70 kHz, and i need to lower them to about 20 kHz (human hearing range 20 Hz - 20 kHz), about 4 times lower.
If i simply change the sampling rate to 250000*4, i get the right frequency but the file become shorter/faster.

I have tried a few things:

Delete every 4th bin of the original vector and multiplying sampling rate by 4. it did the opposite of what i intended and also changed the length.
Duplicate every bin the the audio 4 times and keep the sampling rate, it actually lowered the frequency and the length did not change, but the frequencies were mirrored it 4 times
A method i found online,which i don't understand, this method gave the best result: lower frequency and same file length but it changed the structure of the sound file (see link to spectrogram)

here is the code sample of the last method:

function [signal]=lower_the_freq(original)
% original = rawClip.values';
    binlen=length(original);
    sfq=250000;
    duration=binlen/sfq;

% plot original signal
     subplot(211);

 fourierTransform = fft(original);
 spectrogram(original,1024,512,1024,sfq,'yaxis');    %run the spectrogram
 title('Original signal')

% downsample spectrum by a factor of 2
     n = 2; % downsampling factor
     newSpectrum = fourierTransform(1:n:end);

% zero-pad the positive and negative ends of the spectrum
    pad = floor(length(fourierTransform)/4);
    fourierTransform = [zeros(1,pad) fftshift(newSpectrum) zeros(1,pad)];

% inverse transform
     signal = ifft(length(original)*fftshift(fourierTransform),'symmetric');

% plot the downshifted signal
     subplot(212)
     spectrogram(signal,1024,512,1024,sfq,'yaxis');    %run the spectrogram
     title('Shifted signal')

image


texts:
I'll be processing stereo music signal. What I need is to first split the frequency range into 2 bands at 50-100 Hz. In the higher band, I need to delay one of the channels by a certain number of milliseconds. The lower band has to be delayed as a whole (left+right channel) by a different value. Is this doable? Is this doable in real-time (think VST plugin for an audio player)? If it's possible, can you please give me pointers (algorithm-wise) on how to implement it?

P. S. I need to preserve the sound quality. If such a filter would ruin it, then my idea goes to the trash can.


texts: Spreading symbols with sequence before transmission

texts:
I  have been learning IEEE 802.11ad and the following processing algorithm

I was wondering what is the need for spreading with Golay sequence. I am guessing it has to with something like spread spectrum in CDMA applications.

Does anyone know why and whats the purpose of using such a block after modulation block?

Thanks


texts: Establishing Phrase Frequency in Audio Files

texts:
Is there any Application, Programming Library or Theory that could be used to identify human language phrases in a set of audio files, so as to calculate the frequency of occurrence of each such phrase. For example, how often does the phrase "Thank You" occur in these files?

One solution might be to convert the audio to an International Phonetic Alphabet form but I cant find such a tool.

I'm basically trying to find patterns in speech and the frequency at which they occur. Of course the pattern for "Thank You" is slightly different in different dialects, but I've no problem in having different matches in such cases.


texts: MFCC extraction for speaker verification

texts:
Regarding extraction of MFCC from speech signal...

I have checked for the detailed procedure of extraction of MFCC from http://practicalcryptography.com/miscellaneous/machine-learning/guide-mel-frequency-cepstral-coefficients-mfccs/

For each frame they have extracted the features..
If the signal consists of 100 frames..
then each frame will give 39 values as such..
there by we will get a 39 X 100 matrix for a single speech signal

My question is how to get 39 MFCC coefficients for the whole single signal?
People are saying that they have used 39 MFCCs for a signal, How?


texts: Can you explain me in a easy way what the use of Cyclic Prefix in OFDM is?

texts:
I cannot understand why cyclic prefix is used. I know that it mitigates the effect of ICI but I can't see how. Can you explain please? Formulas and mathematics is well accepted. Thanks.

I read this Cyclic prefix vs. guard interval. Why use one over the other? but it doesn't explain me why. I still do not understand


texts: Getting stuck into an open source image processing/computer vision project

texts:
I have a degree in Electrical Engineering (focusing on VLSI). However, I am more mathematically minded, and only recently discovered that image processing seems to contain lots of interesting maths. I am looking to break into this field in a about a year or so in spite of not having any formal experience in this area.

I was wondering: would contributing to code development on one of the open source projects help my cause? Would this enable me to overcome my lack of formal education or experience in this field?

I would also appreciate it if someone could point me to good open source resources. I know of GitHub and Nasa's Vision workbench. However, I have not contributed to any open source projects and I don't know what to expect.


texts: Detect non rectangular bounding box for objects

texts:
I have a short video where a TV is shown. I need to detect the TV there (let's assume it's all pure black or it's covered with a green screen), and overlay a video on top of it on the same area. But the problem is when the angle is not straight, it is oriented and the bounding boz is not a rectangle. So I guess I need to detect the angle and its scale to overlay the content correctly and nicely.
Is there any ways to do it with AI or computer vision? We can detect the TV with object detection, but I don't know how to compute the tilt and angle for best fit when it's not a retangle.


texts: Implementing a beamforing algorithm for radar

texts:
I want to implement Real-time Radar Beamformer algorithm in C language, and the Conventional beamforming algorithm is as follow :
for each (j:=doppler pulse x range)
   for each (i:=beam)
      sum :=0
      for each (k :=channel)
         sum += c[i][k] * x[k][j]
         y[i][j]:=sum
   end
end

where the x=(x1,x2,..,xns) is the received signal values from n sensors, the y=(y1,y2,..,ynb) the formed beams, and the c is a matrix that composed of the beamforming complex coefficients. So the question that I have is the x will be a vector or a matrix and why c is a matrix ?


texts: 2-D Wavelet output changing when using 'figure' command in Matlab

texts:
I am having trouble understanding what is causing the problem to my Wavelet output when i am trying to display the image with the figure command.

I want to display both the Original image and its Wavelets decomposition,
but when I'm using this code the original image is not displaying in a separate figure.

Here is my code,

im = imread('cameraman.tif');
imshow(im),title('Original Image')

[cA1,cH1,cV1,cD1] = dwt2(im,'bior3.7');
A1 = upcoef2('a',cA1,'bior3.7',1);
H1 = upcoef2('h',cH1,'bior3.7',1);
V1 = upcoef2('v',cV1,'bior3.7',1);
D1 = upcoef2('d',cD1,'bior3.7',1);

% To display the results of the level 1 decomposition, type:

subplot(2,2,1); image(wcodemat(A1,192));
title('Approximation A1')

subplot(2,2,2); image(wcodemat(H1,192));
title('Horizontal Detail H1')

subplot(2,2,3); image(wcodemat(V1,192));
title('Vertical Detail V1')

subplot(2,2,4); image(wcodemat(D1,192));
title('Diagonal Detail D1')

and the output of this code is just one figure with four plots (Digonal, Horizontal, Vertical details and Approximations)

but when I'm using figure command to display the original image in a separate figure,

 im = imread('cameraman.tif');
imshow(im),title('Original Image')
figure

it messes up the wavelet output,

Why does it happens that my output is changed to blue and red color ?


texts: What is best practice to remove static noise?

texts:
What is best practice/algorithms  to remove static noise from speech ? static noise aka Stationary Signal/static noise due to

1.humming sound due to of telephonic line/channel,
2.moving Vehicle,or due to
3.constant noise of fan or other nearby instrument .


texts: Possibilties of chebyshev polynomial waveshaping

texts:
When talking about harmonic distortion, or more specifically waveshaping, we say the order of distortion, can be solved from the equation:

$x^n$

Where n is the order of distortion. To my knowledge the complete result of any waveshaping function can be had by summing all of the different orders of harmonic distortion, multiplied with an appropriate index $a_n$.

$$\sum_{n=1}^{\infty} a_nx^n = WaveshaperFunction$$

Be aware that by feeding the equation with $x=sin(\omega t)$, the highest frequency any $x^n$ can produce is $sin(n*\omega t) $. However, higher orders can still produce multiple different frequencies for a sinusoidal input.

Chebyshev polynomials addresses the issue, removing any other frequencies limiting the harmonic distortion to only $sin(n*\omega t)$, producing an orthogonal result. The polynomials are defined as:

$T_0(x) = 1$

$T_1(x) = 2x-1$

$T_{n+1}(x)=2xT_n-T_{n-1}$

Now, since the result of putting any $sin(\omega t)$ through the chebyshev waveshaper is a sine with a frequency multiple of N, it is easy to say that any harmonic signal with the original signal's phase (or inverted phase) can be produced by summing over all the chebyshev polynomials. That is to say, summing over all the polynomials, any waveshaping function can be generated.

$$\sum_{n=1}^{\infty} b_nT_n(x) = WaveshaperFunction$$

Where $b_n$ represents the amplitude of nth generated harmonic, given a sine input. $T_n$ represents the cheyshev polynomial producing nth harmonic given sine input. I think there could be many advantages for thinking about time-invariant harmonic distortion or intermodulation distortion through this framework rather than through the before-mentioned orders.

As for the question, did I make a mistake somewhere? Particularly I am not so certain that any function can be generated through summing the chebyshev polynomials. Chebyshev polynomials seem to not care about the amplitude of the input sine wave (apart from some DC offset differences), yet I can easily generate a waveshaper function that dose so (any function that is linear near 0 and then changes). Re-framed, can the polynomials generate any possible function?


texts: Can we do cyclic prefix with zeros in OFDM rather than with samples from end of the ofdm symbol?

texts:
considering an OFDM symbol of x(0)x(1)x(2)x(3)x(4)....x(n-2)x(n-1)x(n).
To perform cyclic prefix we prefix some samples from end of this OFDM symbol to the beginning of the OFDM symbol,  like

x(n-3)x(n-2)x(n-1)x(n)x(0)x(1)x(2)x(3)........ x(n-4)x(n-3)x(n-2)x(n-1)x(n).

Instead of taking x(n-3)x(n-2)x(n-1)x(n), if we prefix with zeros what will happen?


texts: Are there any pictures of a sample accurate spectrogram on the web or in books?

texts:
I programmed a sample accurate spectrogram some time ago and i still have never seen a picture of another one. It's interesting maths, i was told it's trivial, so why has no one else posted research on sample accurate sgrams on the internet? Can you find one. Here is a pic of my version:

this is the same program without the zoom and in 2d Mode:


texts:
I am a bit confused as to how to go about this.

I have a surface [x,y,z] given by:

[x, y] = meshgrid(-10:10,-10:10);
z = x .* exp(-x.^2 - y.^2);
surf(x, y, z);

I want to use $I=\rho n^TS$ where $\rho=0.1$, $S$ is the sources vector and $n$ is the normal to the surface, to reconstruct the image. Say I want to do this for 10 arbitrary $S$ vectors and reconstruct 10 images. how do I do this?

I tried extracting $n$ by [nx, ny, nz] =surfnorm(x,y,z); but I am confused here as to how to represent the $S$ vector. is it supposed to $2\times 21$? for one vector? Is there an easier way I am missing?


texts: energy normalization across different scales in case of discrete wavelet transform

texts: In case of continuous wavelet transform (CWT), the wavelets are generated from the mother wavelet by scaling and translation. To achieve energy normalization and to ensure that all wavelets have the same energy regardless of their scales, each wavelet is divided by the square root of the the scale S. In case of discrete wavelet transform, as shown in the attached figure, the energy of the wavelet having 4B bandwidth at the first scale is twice the energy of the wavelet having a bandwidth 2B at the second scale, and so on. Thus, the energy is not equal at different scales. How energy is normalized at different scales in case of DWT?


texts:
I need gaussian pulse fft. I set the sigma of gaussian pulse 1e-10 and define the gaussian function in time and use the fft. as you know the amplitude of should be unit in frequency domain. but my result shows that the amplitude of the pulse in frequency domain equal to 1/end_time. for example, if i define the time as t=0:1/fs:1e-07 then the amplitude of X(f) is 1e7 and I changed the end_time then the amplitude respectively change.

this is my code :

clear all;
clc;
fs=1e10; %sampling frequency
sigma=1e-10;

end_time=1e-7;

t=0:1/fs:1e-7; %time base

t_s=1e-8;
variance=sigma^2;
x=1/(sqrt(2*pi*variance))*(exp(-(t-t_s).^2/(2*variance)));
subplot(2,1,1)
plot(t,x,'b');
title(['Gaussian Pulse \sigma=', num2str(sigma),'s']);
xlabel('Time(s)');
ylabel('Amplitude');

L=length(x);

X =fftshift(fft(x,L));
Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
f = fs*(0:L-1)/L; %Frequency Vector

subplot(2,1,2)
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([0 1e10])


texts: Measuring length of DNA fibers from an image of single molecules

texts:
I am a biologist with very little experience with image processing but have sufficient knowledge of MATLAB and have the image processing toolbox. Ideally I am looking for a MATLAB based solution, but an approach outlining how to go about it would also be helpful.

Update (28 Nov 2011)
It appears that there are certain problems (such as overlaps in signal and definition of color) when using composite images (which is what I presented in the initial question).  I am attaching separate images from the 2 channels: green  and red  (the turquoise regions in the composite image can be ignored), and the coposite image . The red channel is bad for 2 reasons: 1. It has poor contrast due to higher background, 2. Since Red seems to bleed into the green at the background level.

A feature is defined as a region on the composite image that has
Green-Red-turquoise-Red-Green or equivalently the 2 adjacent linear segments on the green and the red that are colinear and contagious.

I am hoping that looking at the images from the two separate channels makes identification of the features easier.

I have the following suggestions for the algorithm:

First identify co-linear green segments (and determine lengths of the green segments)
Determine if there are adjacent contagious and colinear segments facing toward each other (i.e. green->red->  <-red <-green) in the red channel. If yes define the red segment length from the point where green segments end  (because they will overlap with the green segments) till the point on the red segment that is closest to the other red segment of the feature. (i.e. one of the ends of the red segment is set to the end of the overlapping green segment).

Many thanks!

Background:

My question relates to extracting feature from an image:

The original image (tif) is located here:

  Image example 1 (dropbox)

This image is a composite of 3 channels (in tif format): red, green and turquoise. The turquoise colored fibers simply marks all the DNA we have on the coverslip. The feature of interest is the Green-Red--turquoise--Red-Green feature on the single DNA strand that is the middle of the image.

Red is generally the noisiest. This example is good because the contrast is good. However, sometimes the images are not so nice and there is hue throughout the image, so hard-coding a specific RGB value for the green and red color might not work for all images. Also, note that the fibers are not necessarily horizontal, they might be rotated (but never vertical).

Please see this image for an example:

The original image (tif) is located here:

  Image example 2 (dropbox)

Also, sometimes a single image has many such features and sometimes there are multiple features on the same DNA strand. Finally sometimes there might be only partial features (i.e. isolated green or isolated red or isolated green-red segments, but unpaired).

Question:

I would be grateful if someone can help me obtain the lengths of the individual segments of green and the red segments i.e. since the feature of interest is Green-Red--turquoise--Red-Green, each feature would have an array of 5 values (length of the first green segment, length of the first red segment, length of the turquoise segment, length of the second Red segment and the length of the second Green segment).


texts:
Is it true or fale:

Given LTI, $H(z)$, with impulse response $h[n]$ that is not zero for all $n$. The following system is not LTI:

$$x[n]\longrightarrow\boxed{\downarrow2}\longrightarrow H(z)\longrightarrow\boxed{\uparrow 2}\longrightarrow y[n]$$

TRUE

Why? If I take delta to be the input, and then take another input: delta shifted by one...how I see that it will support the statement?


texts: Blind source separation with one source rotated

texts:
I have two images:

$$
I_1 = w_{11}A + w_{12}B, \\
I_2 = w_{21}A + w_{22} \bar{B}.
$$

$A$ and  $B$ are unknown. $\bar{B}$ is ${B}$ rotated by 180 degree. For both images, $A$ has higher signal-to-noise ratio. What would be an effective way to extract $A$ out of these two images?


texts: How to obtain the steady state response of the sampled periodic signal data?

texts:
Below the green plot Vin is 5Hz sine voltage signal from a function generator output sampled at 12kHz sampling rate:

The blue plot Vout is the processed/filtered data by using the following code in Python:

#Transfer function
fc = 1 # roll off freq
wc = 2*numpy.pi*fc
sys = signal.lti([1], [1/wc, 1])  # transfer function

#Output in time (Response of the filter to the input)
tout, v_out, x = signal.lsim(sys, v_in, time)  # output

And here is Vout in a more detailed view:

Basically I define a low pass filter transfer function  using:
 sys = signal.lti([ 1 ], [1/wc, 1]).

Then the Vout is obtained from tout, v_out, x = signal.lsim(sys, v_in, time)

So Vout is the filtered response of Vin.

If you compare the Vin and Vout in above plots, Vout looks like it is processed like in Laplace domain, that is to say Vout looks like it is both the transient and steady state response of Vin.

(This problem seen more obvious for a DC input; the output Vout becomes like a transient response of an RC filter)

Having the sampled data like in green plot Vin, how can we obtain a steady sate response of such filter? Or at least close to steady state.


texts:
I am measuring wind speed using a number of different anemometers. These anemometers are placed on a structure which exhibits sinusoidal motion. Thus, the signal of the anemometers contains both the wind speed as well as the sine wave to which they were exposed.

I wish to accurately determine the amplitude of the sine wave from the anemometer signal. Thus, I decided to use a flat-top window when calculating the FFT. However, since the signal is 'noisy' the resulting amplitude of the sine wave is not correct. Sometimes, the amplitude given on the spectrum is greater than the actual amplitude of the sine wave!!

I am considering using Welch's overlapped method (pwelch in MATLAB). Can someone please tell me whether such a method is suitable for my application? I have read that this method should only be used with stationary signals? I am confused whether my signal is considered as stationary. Can someone please help me? Thanks.


texts: How to perform model fitting for system identification

texts: I am having a really hard time in understanding how to formulate a model say linear AR model to represent a communication channel or maybe any motion. I have the experimental data representing the kinematics of robot movement. In the case of a communication channel, it will have a medium, transmitter and receiver. So, will the medium be represented by AR model or do we represent the transmitter and receiver by these models? How to decide the order of the system to be chosen initially. I am not from the area of signal processing but need to use this for a part of my research work. I have gone through literature review and the background of Kalman filter, Least square method, recursive least square for estimation purpose. But I cannot find a good example which shows how the parameters are estimated from the time series or motion data or speech modeling. Can somebody explain or point out the starting point and how to go about it? I know this is a very broad topic but an explaination with a small example will be really an eye opener for a beginner.


texts: Scale and Rotation invariant Template Matching

texts:
I'm looking for a method for scale and rotation invariant Template matching. I already tried some, but they didn't work so good for my examples or took for ever to execute
.
SIFT and SURF Feature detection failed totally. I also tried to implement a Log-Polar Template Matching function, but I never finished (didn't know exactly how to).

In these Articles (the first is in german)

http://cvpr.uni-muenster.de/teaching/ss08/seminarSS08/downloads/Wentker-Vortrag.pdf

http://www.jprr.org/index.php/jprr/article/viewFile/355/148

I read about that method. Mapping the polar coordinates worked, but I don't know if it's right. The Images look like this.

source_log_polar.png http://www.shareimages.com/images/pics/0/0/3/62394-pZSfl5WenZysnpyVnKg-source_log_polar.png

And after matching these 2 images with OpenCV's Template Matching function i got that result

Now I don't how to go on.

My templates are always simple symbols in building blueprints and the blueprints itself. The symbols can differ in size and orientation.

For example my simple blueprint:

And my template

In this example there is only one template, but in the blueprints it should find all occurrences, even the ones with sizes and/or orientations.

Does anyone have an approach how I could solve this?

Edit:

An addition to Andrey's approach. The distance capturing algorithm for a radial profile. (Using EmguCV)

private float[] getRadialProfile( Image<Gray, byte> image, Point center, int resolution )
 {

 var roi = image.ROI;

 if ( !roi.Contains( center ) )
  {
   return null;
  }

 var steps = resolution;
 var degreeSteps = 360 / (double)resolution;
 var data = image.Data;
 var peak = 0.0f;
 var bottom = double.MaxValue;
 var bottomIndex = 0;
 var width = roi.Width;
 var height = roi.Height;
 var minX = roi.X;
 var minY = roi.Y;

 float[] distances = new float[resolution];
 for ( var i = 0; i < steps; i++ )
  {
   var degree = i * degreeSteps;
   var radial = degree * Math.PI / 180.0;
   var dy = Math.Sin( radial );
   var dx = Math.Cos( radial );

   var x = (double)center.X;
   var y = (double)center.Y;

   while ( true )
    {
    x += dx;
    y += dy;
    if ( x >= minX + width || y >= minY + height || x <= minX || y <= minY )
     {
      x = -1;
      y = -1;
      break;
     }
    var pixel = data[(int)y, (int)x, 0];
    if ( pixel == 0 )
     {
      break;
     }
    }

    float distance = 0.0f;
    if ( x != -1 && y != -1 )
    {
      distance = (float)Math.Sqrt( Math.Pow( (center.X - x), 2 ) + Math.Pow( (center.Y - y), 2 ) );
    }

    distances[i] = distance;
    if ( distance > peak )
    {
      peak = distance;
    }
    if ( distance < bottom )
    {
      bottom = distance;
      bottomIndex = i;
    }
   }

    // Scale invariance. Divide by peak
   for ( var i = 0; i < distances.Length; i++ )
   {
     distances[i] /= peak;
   }

    // rotation invariance, shift to lowest value
   for ( var i = 0; i < bottomIndex; i++ )
   {
     distances.ShiftLeft(); // Just rotates the array nothing special
   }

   return distances;
}


texts: is root mean square possible to calculate similarity of FFT result?

texts:
I need to compare FFT result from audio files,
I have 2 FFT result from 2 audio files..

FFT 1 (sample)

-0.16314493488504767 - 0.103707391105263i
0.07863935536550609 - 0.30111206509352917i
0.06753127619272284 + 0.02837438103569126i

FFT 2 (sample)

-0.006953384463464084 + 0.0291883094944081i
-0.09447919915711694 - 0.22676541801892i
-0.23493813662827812 - 0.07765408995141115i

Then I need to calculate the deviation of each real-part point (without imaginer part) from the FFT result then make the result absolute, something like this

-0.16314493488504767 - (-0.006953384463464084)
--> then use Math.abs to get absolute value

then sum all the deviation point, and use Root Mean Square formula to follow with to get the difference value...

where N is my numbeofFFTPoint and the following x is the sum value of my deviation result.

But then, someone told me to use Zero crossing rate formula to get more accurate result

Honestly I don't really understand about Zero Crossing Rate,

what have I tried :

I compare the exact same 2 audio files, then get the deviation and sum result --> 0.0  and of course the Root Mean Square Formula result is 0.0 too..
since there's no difference between the files the it would return 0.0 right ??

question

is it possible just to user Root Mean Square ??
Any advice which formula should I use for this case ?

Thanks :) (sorry 4 bad explaining)


texts: Rapid changes in signal correspond to high frequencies - proof

texts:
If I consider a generic aperiodic signal $x(t)$, how can I prove that rapid changes in signal correspond to high frequencies?

Thank you for your time.


texts: Good metric for qualitatively comparing image patches

texts:
I am trying to "match" little square patches in an image. At first glance, it seems reasonable to simply do a Euclidean distance style comparison of two of these arrays to get a "similarity" measure. This works fine in many cases (the "best" patch (lowest value) according to this metric looks very much like the query patch). However, there are many cases in which this produces a very bad match. For example, take these two patch pairs:

Two patches of a brick wall, score 134 (this is the sum of the components of the average absolute pixel difference):

One patch of a brick wall, one patch of grass, score 123!

To a human, "clearly" the grass does not match the brick, but this metric says otherwise. The problem is just in the local statistical variation.

If I use something like a histogram comparison, I completely lose all spatial information - e.g. if a patch is grass on the top and brick on the bottom, it would match exactly to a patch with grass on the bottom and brick on the top (again, another "obviously wrong" match).

Is there a metric which somehow combines both of these ideas into a reasonable value that will evaluate to "similar" for Pair 1 above, but not also be similar for my example of a patch and its vertical mirror?

Any suggestions are appreciated!


texts:
Maybe this a stupid question, but I came along with this point and I would like to confirm that I am not wrong.

I found some similar posts in StackExchange, like:

Expectation and Auto-correlation of an Output from a Matched Filter.

In one of the answer, it is stated "as S(t) is determistic[...]" and then the deterministic signal is taken out of the expectation.

and we know that for instance if I have a constant value, I can take it out of the Expectation operator, e.g.:

Where lambda is a random variable and "a" is a constant. Ok, so, if "a" is a constant, then fair enough, but what happens when "a" is a deterministic signal, but it is not a constant, for instance, what happens if "a" is a squared subcarrier?

Can I still safely say that the expectation will be the product of the deterministic signal multiplied by the expectation of the random variable?

Thanks, guys!


texts:
I'm a reasonably competent web programmer and I would like to know how I should go about learning to program a software that would control a music sampler.
The hardware setup I was given is a sampler
( Example : https://en.wikipedia.org/wiki/Akai_S3000XL or http://www.vintagesynth.com/yamaha/a3000 )
and a PC, but this can be extended.
Somebody should be able to operate the sampler from anywhere by means of a web page that shows a sort of a control panel which should be identical in function and form with the one of the sampler. Then they should save / download the result.
I'm not sure this is the right place to ask, but I don't have any other ideas.
Software wise I can do anything on the webstack (PHP, Javascript, Python, .Net), and also Java, and C++ if need be. I have no experience with Rust, and I won't touch Assembly.
Hardware wise, what kind of hardware setup should I be looking for? This is what I don't understand actually: how to interface with the samplers.
If there's anything like an idiot's guide to control some random machine from a PC, please point me to it.
Thanks for any help, it is much appreciated.


texts: Must OFDM signal be time-division in a multiple-transmitter and SFN system when the signal is free of spread-spectrum?

texts:
The OFDM baseband signal can be demodulated since  each subcarrier is orthogonal to another. But In SFN system, i.e., Single Frequency Network, all OFDM signals from several transmitters are in the same frequency domain, which breaks the orthogonality.

So without some technology like CDMA, the OFDM sigals from each transmitter must be time-divided.

Is that right?


texts: Daubechies Wavelets in mulitresolutional analysis decomposition

texts:
I have an understanding problem with Daubechies wavelets.

When I use a multiresolutional analysis, I want to approximate the given input Signal $f\in L^2(\mathbb{R})$ on the subspaces $V_i$. By calculating the coefficients I split the signal in the detail part and the rest. Now I am wondering, how the Daubechies wavelet is fitting in this picture. The Daubechies wavelet is defined by coefficients, but in my opinion, these are only usable if I try to perform a wavelet decomposition on input data in $l^2(\mathbb{Z})$.

To specify my question, I have no idea how to calculate the first inner products
$$c_i = \langle f, \varphi_{i,j} \rangle, \text{ and } d_i = \langle f, \psi_{i,j} \rangle $$
if $\varphi$ and $\psi$ are Daubechies scaling function and wavelet, since they have no expression as functions in $L^2(\mathbb{R})$.

Thanks in advance

Matthias


texts: Literature on noise removal in OCR applications

texts: Could someone point out the core texts or articles most useful on techniques for the removal of noise from scanned text for OCR applications?


texts:
The deconvolve(signal, divisor) function of scipy

> Deconvolves `divisor` out of `signal`. Returns the quotient and remainder such that `signal = convolve(divisor, quotient) + remainder`

The core codes of deconvolve are

num = atleast_1d(signal)
den = atleast_1d(divisor)
N = len(num)
D = len(den)
input = ones(N - D + 1, float)
input[1:] = 0
quot = lfilter(num, den, input)

It's quite different from any deconvolve methods I ever know, such as the frequency spectrum division.

Can anyone explain those behind the codes?


texts: How to implement a FIR filter with correction of group delay on c#?

texts:
I read respect to digital filters, and implemented both IIR and FIR filters and I found that this type of filters have a called "group delay". Now I know that the IIR have a non-linear group delay and the FIR have a linear group delay.

I have a real-time system that require a high-pass filter computed over a computer (good processor and RAM), allways that pass through the filter I'll have a phase distortion, taking advantage of the good processor: What is the best method to correct the phase and faked the response of the filter with the correct phase?

In the image above I extend my input signal with the same number of coefficients at the right and I try to correct the group delay, filtering the first time, invert the response in the time and pass again by through the filter and revert again the response in the time. The result correct the phase but causes other type the distorsion in the last samples. What is wrong? My knowledge about filters is not so wide unfortunately.

Note: The number of coefficients that I am using is 1000, Fs: 6400 samples/second and cutOff: 10 Hz.

I appreciate so much any help and aclaration. Thanks


texts: How can I distinguish between two similar EEG signals?

texts:
I have two EEG signals that are very similar. The difference is only in amplitudes. However, they are coming from two different cognitive processes. What are some methods, beside FFT, for distinguishing between two very similar signals in terms of statistical information? Please let me know. Thank you

BSS.


texts:
There are at least two ways of achieving sub-sample delay. One is to simply over-sample the signal and then delay it by the wanted sub-sample amount. The over-sampling is equivalent to creating new samples in-between the existing samples utilizing the prior sample values and then using a low-pass filter on the resulting signal. Low-pass filters generate artifacts, resulting at minimum in pre/post-echo for linear phase filters.

The second alternative is to perform a FFT, multiplying the phase response with a linear phase response corresponding to the delay amount. My question is: what are the artifacts involved in this procedure? I would think that they would be the same, but can't find a low-pass filter anywhere.


texts:
I'm using 'awgn(input_signal,snr_db,'measured')' function in MATLAB to add white noise to the input signal.
Here is the simplified script.
Does zero padding the input signal affect the behavior of awgn function?
clc;clear all;close all;

data = randi([0 3], 1, 2^12); % Generate data
qam_mod = qammod(data, 4);  % Modulate data
qam_mod_u = upsample(qam_mod,64); % Upsample

ch_out_interplt = awgn(qam_mod_u, 5, 'measured');
ch_out_no_interplt = awgn(qam_mod, 5, 'measured');

ch_out_downsample = downsample(ch_out_interplt,64); % Downsample

scatterplot(ch_out_downsample); title('with up, down sampling, 5 dB SNR');
scatterplot(ch_out_no_interplt); title('no up, down sampling, 5 dB SNR');

The scatterplots are different with and without zero padding.

Does that mean awgn function behavior depends on the zero padding (length of the signal)?
How does the 'awgn(input_signal,snr_db,'measured')' function calculates the input signal power?
How to make the power calculation independent of the length of the input signal?
Is there any better way to simulate the awgn noise in MATLAB?


texts: Estimating degree of motion in video frames

texts:
I have a video (in webm format), and I want to produce a graph that shows the degree of motion in this video as a function of video time. I'm sure there are many different ways to define "degree of motion", but I envisaged something that could be based on the average difference between pixels in consecutive frames, or something like that (suggestions are welcome). Essentially what I need is a graph whose magnitude is large in cases where there is a very big change in the frames and small when most of the frame stays the same except for a small section.

What tools might I use to achieve this? Thanks in advance.


texts: Time-Frequency analysis of IQ data in MatLab

texts:
I am new to signals processing, and am trying to match the information I'm getting out of MatLab to that taken from a commercial spectral visualization/playback tool. That is, I know what my data should look like in the time and frequency domain, but am thus far unable to correctly compute the latter.
My datasets are IQ recording, for which I have the following additional data: Center frequency, frequency span, and Sample Rate; as well as supplier-provided equations (with constant values) for converting the raw IQ data to complex Voltage values and that Voltage data to real dBm values.
If I want to display a spectrogram wherein the power is accurately represented in dBm, what is the order of operations I need to follow? Thus far, all my attempts have given me wildly inaccurate results (either the power is too low or in excess of the time-domain max amplitude in dBm).


texts: Calculate phase lag between two signals with perturbed frequencies

texts:
This type of question has been asked quite a few times on this forum and others now, but I still haven't found a satisfactory answer to my problem.

Given an input signal: $$x_1(t)=\cos\big(2\pi ft\big)$$
A process derives an output signal: $$x_2(t)=\cos\big(2\pi (f+b)t+\phi\big)$$
My objective is to determine the phase lag $\phi$.

$b$ is an unknown perturbation to the input frequency very close to 0 due to the mechanical nature of the process. For example, a human trying to track this signal.

If the perturbation is exactly zero, I'm able to use cross spectral density to determine the phase lag. However, if the perturbation is even 0.1% of the original frequency, then my answer is totally erroneous.

For example, in Python:

import numpy as np

Fs = 1000
t = np.arange(0,1,1/Fs)
x1 = np.cos(2*np.pi*100*t)
x2 = np.cos(2*np.pi*100.1*t + np.pi/4)
dft1 = np.fft.fft(x1); dft1 = dft1[0:int(n/2)+1]/Fs;
dft2 = np.fft.fft(x2); dft2 = dft2[0:int(n/2)+1]/Fs;
n = len(t)
f = np.multiply(range(0,int(n/2)+1), Fs/n)
csd = np.conj(dft1)*dft2
print(np.rad2deg(np.angle(csd[abs(csd)>0.05])))

Expect 45 but get 62.9 instead.

I've tried Hilbert transform and straight fft. Any idea?


texts: Disagreement between Analytic Convolution and FFT Convolution

texts:
I've been working on a small project, and as a part of this I need a function that can compute the structure function of the input data. In the process of writing this part/testing it I found the results were actually a poor match to what I would expect based on an analytic solution. Upon some further investigation, I've found this mismatch applies even to 2D convolutions and I have been completely stuck trying to debug this.

The test case that I am running is two rectangular functions convolved with one another, my expected result is a 2D triangle function (a pyramid essentially). But the difference between the analytic output and the FFT convolution output is large enough that I am concerned something is wrong. Furthermore, when I compute a direct convolution of the two inputs using scipy.convolve I get yet another output, different from the other two in it's own unique way!

To demonstrate the differences I'm seeing, I've plotted: my FFT Convolution's output, the analytic solution, the difference between them, and the difference between scipy.convolve's output and the analytic solution. My question then boils down to a couple things:

How large of a difference is expected between the FFT convolution and the analytic solution?

I would have expected something small (of order 1e-15 or less) but the peak difference I am seeing is 0.12. It also follows a very particular pattern, the vast majority of the output is close enough (difference of ~1.0e-16).

Is there a key step or concept that I am missing here?

I've been taking the FFT of each input, multiplying them together, and then taking the inverse FFT (with some scaling based on pixel size).

The code that produces the figure can be found in this relatively simple jupyter notebook on my github.


texts:
Prove that the following system is invertible.
$$y(t) = \mathcal{T}\{x(t)\} = \int_{-\infty}^{3t} x(\tau) \,\mathrm d \tau$$

Answer: yes, the system is invertible. I need some hint here, not the full solution.


texts: Poles and zeros representation in s-plane

texts:
We have a filter (or a system) that has two poles at $f_{p1} = f_{p2} = 20$ GHz and one zero at  $f_{z1} = 15$ GHz.  Let's assume the gain is $K = 1$. Then how do you mark these poles and zero in the s-plane?

Is it right that we consider it as three points at the $j \omega$ axis with real part equal to 0? i.e:

$s_{z_1} = 0 + j2 \pi f_{z_1} = j \omega_{z_1}$,

$s_{p_1} = 0 + j2 \pi f_{p_1} = j \omega_{p_1}$, and

$s_{p_2} = 0 + j2 \pi f_{p_2} = j \omega_{p_2}$

So in this case if I want to write the transfer function it would be

$$ H(s) = \frac{s - j \omega_{z_1}}{(s - j \omega_{p_1})(s - j \omega_{p_2})}$$

Then I can derive $H(j\omega)$ by substituting $s$ with $j \omega$

$$ H(j\omega) = \frac{j\omega - j \omega_{z_1}}{(j\omega - j \omega_{p_1})(j\omega - j \omega_{p_2})} = \frac{\omega - \omega_{z_1}}{(\omega -  \omega_{p_1})(\omega - \omega_{p_2})}$$

Is that right so far?

Now if I want to convert the zero location into discrete time domain then which of the following should I use:

a) $\frac{1+\frac{\omega_{z_1}}{f_s}}{1-\frac{\omega_{z_1}}{f_s}}$

or

b) $\frac{1+\frac{s_{z_1}}{f_s}}{1-\frac{s_{z_1}}{f_s}}$ = $\frac{1+\frac{j\omega_{z_1}}{f_s}}{1-\frac{j\omega_{z_1}}{f_s}}$


texts: Do I need to bandpass filter the signal before applying the Zoom FFT algorithm?

texts:
Let's say I have an acoustic signal that I want to make a spectrogram of filled with broadband noise and tonals, sampled at oh, $44100\textrm{ kHz}$, and I want to zoom in on the band between $2600\textrm{ Hz}$ and $3000\textrm{ Hz}$ (or any band of frequencies, really).

The zoom FFT algorithm says to accomplish this, one simply has to frequency translate, then lowpass filter, then decimate and FFT the resulting lower frequency signal.

However, most of the diagrams I see explaining zoom FFT only show discreet bands with no noise in between, which begs the question:

Would I need a bandpass filter to prevent aliasing during the frequency translation step of the Zoom-FFT algorithm if I have a noisy signal containing energy at all frequencies?

I'm thinking yes, since there's a chance that during the frequency translation some unwanted noise would get shifted into my new band of interest due to aliasing, but I'm not sure.

The processing chain I envision is as follows:
$$\boxed{\text{BPF}}{\rightarrow}\boxed{\text{Frequency Xlate}}{\rightarrow}\boxed{\text{Bandlimited downsample}}{\rightarrow}\boxed{\text{FFT}}{\rightarrow}\boxed{\text{Display processing}}$$

Would this be correct?


texts: How do I get texture feature images from GLCM texture analysis?

texts:
In this paper under section 2. "Proposed Methodology", the author creates a GLCM from an image and then extracts texture features resulting in a new image for each feature.  I've seen other authors do this in other papers as well.  The texture features used in the image are the common Haralick features detailed on page 619 here.

How did the author extract a texture feature image from those features?

The author has a unique image for contrast, an image for entropy, an image for correlation, etc.  When I look at the Haralick feature equations the result is always a number not another image.
How is a unique image created for each texture feature?

EDIT:

For GLCM definition as suggested in comments, a GLCM stands for the Gray-level co-occurence matrix.  It's a common technique used for texture analysis.  It is defined over an image to be the distribution of co-occurring pixel values at a given offset.


texts: Poisson summation formula and periodic summation of Fourier transforms

texts:
One of the forms of the Poisson summation formula is

$$
\sum_{n=-\infty}^{\infty} T\cdot x(nT)\ e^{-i 2\pi f T n}\;
{=} \;
 \sum_{k=-\infty}^{\infty} X\left(f - k/T\right),$$

where $x(nT)$ are samples of a continuous function $x(t)$, and $X(f)$ the fourier transform of $x(t)$. The RHS is a periodization of $X(f)$.

What I couldn't understand is the following, we can define two functions $x_{1}(t)$ and $x_{2}(t)$ that are equal at the sample points, $x_{1}(nT)=x_{2}(nT)$. We  have
 $$
\sum_{n=-\infty}^{\infty} T\cdot x_{1}(nT)\ e^{-i 2\pi f T n}\;
{=} \;
 \sum_{k=-\infty}^{\infty} X_{1}\left(f - k/T\right),$$
and
 $$
\sum_{n=-\infty}^{\infty} T\cdot x_{2}(nT)\ e^{-i 2\pi f T n}\;
{=} \;
 \sum_{k=-\infty}^{\infty} X_{2}\left(f - k/T\right),$$

where $X_{1}(f)$ and $X_{2}(f)$ are the Fourier transforms of $x_{1}(t)$ and $x_{2}(t)$ respectively.

But since $x_{1}(nT)=x_{2}(nT)$, we then also have $$
\sum_{n=-\infty}^{\infty} T\cdot x_{1}(nT)\ e^{-i 2\pi f T n}=\sum_{n=-\infty}^{\infty} T\cdot x_{2}(nT)\ e^{-i 2\pi f T n},$$

Which means $$\boxed{\sum_{k=-\infty}^{\infty} X_{1}\left(f - k/T\right){=} \;
 \sum_{k=-\infty}^{\infty} X_{2}\left(f - k/T\right)}$$

But we can choose the functions $x_{1}(t)$ and $x_{2}(t)$, even though they are equal at the sample points, in such a way that their transforms $X_{1}(f)$ and $X_{2}(f)$ are quite different from each other, and, therefore, so will their periodizations. And so why should this last equation be true in general?


texts: Find the impulse response from the frequency response

texts:
So I'm having a problem here which gives me the frequency response and asks for the impulse response:

$H(\Omega ) = e^{-j\frac{\pi }{2}}$ for $\Omega>0 $ and
$H(\Omega ) = e^{j\frac{\pi }{2}}$ for $\Omega<0$

I need to find the impulse response. My first thought was to (obviously) apply Inverse Fourier Transform but I'm thinking that since it's a branch function maybe that's not the correct approach. Also, I see that the two forms look a lot like each other so maybe I can somehow combine them into one and THEN apply the Transform. Any thoughts and guidance?


texts: How would I fit a distribution to this image noise?

texts:
I have collected some noise data from a dimly lit CMOS image sensor.  The distribution of pixel values is tallied below:-

I'd like to be able to simulate this sensor noise.  How would I fit a statistical distribution to this bell curvish graph? And which distribution? Bear in mind that the distribution is discrete.  From the raw samples, I get mean = 18 and standard deviation = 7 (both approximations).


texts: GNU Radio Companion - Constellation Modulator - what is the reason behind the interesting Time Sink waveforms obtained?

texts:
For my setup of a Constellation Modulator block and QT GUI Time Sink block, with the symbol map for the Constellation Object being [0, 1, 2, 3], and the Constellation Points for the Constellation Object being [-1-1j, 1-1j, 1+1j, -1+1j] -------- the GUI Time Sink block shows waveform plots that aren't what I expect to see. I can see signal values being -0.707 (which might be ok). But other observed 'signals' show interesting waveforms - that I don't expect to see.
I used GNU Radio Companion 3.8 to run four simulations based on the same project model. The first simulation is based on the Vector Source having a '0' for the vector. Then, I have another run for Vector set to '1', and then 2', and then '3'. That covers the four 'vectors'.
Each screenshot shows Signal 1 and Signal 2 plots. Do I need to fix up my setup for my simulation in order to get expected waveforms - rather than those unusual looking waveforms? Thanks in advance for recommendations/suggestions!


texts: There appears to be a contradiction with the recurrence solutiom

texts:
Consider the equation :
$$
y[n]-\frac{1}{2}y[n-1]=x[n]
$$
where $x[n]:=\left(\frac{1}{3}\right)^{n}u[n]$. First, I founded the homogenous solution $y_{h}:=A\frac{1}{2^{n}}$
$$
y_{h}[n]=\frac{1}{2^{n}}y[0]
$$
where I insist $A=y[0]$. Furthermore, The particular solution is given of the form $y_{p}[n]:=B\left(\frac{1}{3}\right)^{n}$ I shall determine $B$ as follow :
$$
B\left(\frac{1}{3}\right)^{n}u[n]-\frac{1}{2}B\left(\frac{1}{3}\right)^{n-1}u[n-1]=\left(\frac{1}{3}\right)^{n}u[n]
$$
\begin{align*}
\implies B&=\frac{\displaystyle\left(\frac{1}{3}\right)^{n}u[n]}{\displaystyle\left(\frac{1}{3}\right)^{n}u[n]-\frac{1}{2}\displaystyle\left(\frac{1}{3}\right)^{n-1}u[n-1]}\\ \\
&=\frac{\displaystyle\left(\frac{1}{3}\right)^{n}}{\displaystyle\left(\frac{1}{3}\right)^{n}-\displaystyle\frac{1}{2}\left(\frac{1}{3}\right)^{n-1}}\\ \\
&=-2
\end{align*}
So now that we have $B=-2$, We are given the initial rest condition so what I did was sub $n=0$ in difference equation
$$
y[0]-\frac{1}{2}y[-1]=x[0]
$$
$$
\implies y[0]=1
$$
$$
**\text{So this means $A=y[0]=1$}**
$$
However, the major issue is that :
$$
y[0]-\frac{1}{2}y[-1]=x[0]
$$
$$
\implies y[0]=1
$$
$$
\implies y[n]:= y_{h}[n]+y_{p}[n] \implies y[0]=y_{h}[0]+y_{p}[0]=1
$$
$$
\implies A\left(\frac{1}{2}\right)^{0} -2\left(\frac{1}{3}\right)^{0}=1
$$
$$
\implies A=3
$$
So I have $A=1$ and $A=3$, an obvious contradiction and I hope someone can help me


texts: Estimating Convolution Kernel from Input and Output Images

texts:
Given an original image and a convolved version of it, I need to calculate the convolution kernel. For example, given:

and

I would like to calculate the convolution kernel that generated the second image from the first (a gaussian in this case).

The images will be similar to the example, but I can change the configuration slightly to make the process easier (change the shape, have a number of shapes etc.)


texts: What would an algorithm that can target different devices be called?

texts:
I'm developing an audio filtering algorithm that can be run in different scenarios, e.g as embedded firmware on DSP Codec/Core, as a software library, or as a cloud service. Effectively, it can target solutions for hardware, software, and services - but I realize there are semantics and overlap here.
I'm struggling to identify a category to describe this breadth of approach.
An analogous question, by way of example, would be what would I call a software library that can run on Windows, Linux, and macOS or Intel, Arm, etc. The answer would be a cross-platform software.
Can anyone help?
Edit 1:
following on from a comment below, it is useful to focus on the 'portable' nature and it being loosely-coupled to any architecture, device or platform.


texts: I can't understand why there is a phase shift in phase modulation of a sawtooth wave

texts:
This is regarding question 5.1-2 (see picture attached) from Modern digital and analog communications, Lathi & Zhi Ding (2010).

I do not understand what the behavior should be at the discontinuity (sudden jump from 1 to -1). I ran the following code on MATLAB (using simpler values):
a = 0.25;
t = -a:0.0000001:a; % time
wc = 2*pi*1e2; % carrier frequency
m = sawtooth(2*pi*5*t,0.9); % message signal
phi_am = cos(wc*t+(pi/2)*m); % phase modulated signal
plot(t,phi_am); ylim([-2 2]); hold on;
plot(t,m); hold off;

and it showed me a phase shift by $\pi$, and that is what the solution is. But I do not understand how.
Also, if some light on why must $k_p < \pi$ in part (b) then I'd appreciate it.


texts: I am doing some project on medical image processing and I need some uncompressed medical images especially magnetic resonance angiography, vessel and so on. I have searched a lot but most of the available images are in JPEG format which is not a suitable format due to its lossy compression. Do you know any free dataset which have uncompressed medical images especially those mentioned?


texts: If I line in the sound from a video game console, can I detect when a sound is played?

texts: The question is as the question is. It won't be this isolated sound though, it'll be a sound on top of background music. I think I can get a hold of the sound file itself from the data, but at worst I'd just make my own recording and ask for a less accurate match - if variation is possible in whatever method may be proposed.


texts: Why use transfer functions than differential equations?

texts:
I have some simple questions regarding DSP and Control Systems, but found no simple answers. I just need simple and easy to understand answers, not complex answers with huge maths.

Why we use Transfer Functions, when we can get a system's output by just solving it's differential equation?
What we finally get after solving a transfer function (In Laplace Transform)?


texts: Why does increasing duty cycle decrease the separation between spectral components?

texts: I know that the duty cycle measures the width of the signal in percentages. But why does it decrease the signal width with decreasing duty cycle?


texts: Is there an a method to fit a wave created from two wave?

texts:
I need to get the frequency and amplitude for a wave that consists of multiple function.
for example, if I have a sine curve (created from two sine waves), How can I extract the parameters for this complex wave.
I used the curve fitting code provided in this answer, as illustrated below. A simple sinusoidal function is used in the provided code. so I obtain a fitted curve with a single amplitude and frequency.
Data example:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.arange(0, 50, 0.01)
x2 = np.arange(0, 100, 0.02)
x3 = np.arange(0, 150, 0.03)
sin1 = np.sin(x)
sin2 = np.sin(x2)
sin3= np.sin(x3/2)

sin4 = sin1 + sin2+sin3
plt.plot(x, sin4)
plt.show()

import numpy, scipy.optimize

def fit_sin(tt, yy):
    '''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''
    tt = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

    def sinfunc(t, A, w, p, c):  return A * numpy.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*numpy.pi)
    fitfunc = lambda t: A * numpy.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)}

yy = sin4
tt = x
res = fit_sin(tt, yy)
print(str(i), "Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res )
fit_values=res["fitfunc"](tt)
Frequenc_fit= res['freq']
print(i, Frequenc_fit)
Frequenc_fit=Frequenc_fit
Amp_fit=res['amp']
Omega_fit=res['omega']
Phase_fit=res['phase']
Offset_fit=res['offset']
maxcov_fit=res['maxcov']
plt.plot(tt, yy, "-k", label="y", linewidth=2)
plt.plot(tt,fit_values, "r-", label="y fit curve", linewidth=2)
plt.legend(loc="best")
plt.show()

I got a fitted sine curve with a single frequency and amplitude as follows:

Is there a method to obtain fitted curve matches with the original one? or open software can be used to get the parameters for the best-fitted wave from the original one?
PLEASE, my understanding of FFT is VERY limited! So PLEASE forgive me for my lack of knowledge.


texts:
I have below doubts, so confusing! As I don't want to assume from what I read, I am asking for help here!

Are Second Order Sections another name for biQuads ?

If I have 2 single pole transfer functions, then cascading them equals multiplying the polynomial; such as below?

$$H_1 =\frac{b_{10}+b_{11}Z^{-1}}{1-a_{10} Z^{-1}}$$ and $$H_2 =\frac{b_{20}+b_{21}Z^{-1}}{1-a_{20} Z^{-1}}$$

Then cascaded version is a 2nd order polynomial such as:

$$H_c =\frac{b_{0}+b_{1}Z^{-1}+b_{2}Z^{-2}}{1-a_{0} Z^{-1}-a_{1} Z^{-2}}$$ where the poles and zeros can be obtained by taking root of this polynomial, which is not the two individual poles or zeros exist before cascading! Am I right? Or taking square root will give the same poles or zeros of the individual ?

What is the mathematical steps behind the MATLAB commands tf2zp and zp2sos? How to do the same manually [just like we handcode calculation of the FFT using the butterfly tree diagram]?

This shows the Direct form II Transposed structure and it's implementation equation. Now taking values of $a_0$, $a_1$, $b_0$,$b_1$,$b_2$ from $H_c$ and using them in the implementation equation is right?


texts:
My goal is to obtain a plot with the spatial frequencies of an image - kind of like doing a fourier transformation on it. I don't care about the position on the image of features with the frequency f (for instance); I'd just like to have a graphic which tells me how much of every frequency I have (the amplitude for a frequency band could be represented by the sum of contrasts with that frequency).

I am trying to do this via the numpy.fft.fft2 function.

Here is a link to a minimal example portraying my use case.

As it turns out I only get distinctly larger values for frequencies[:30,:30], and of these the absolute highest value is frequencies[0,0]. How can I interpret this?

What exactly does the amplitude of each value stand for?
What does it mean that my highest value is in frequency[0,0] What is a 0 Hz frequency?
Can I bin the values somehow so that my frequency spectrum is orientation agnostic?


texts:
I was hoping that you can help me understand how to find poles and zeros in this transfer function.

$$
H(s)=\frac{1}{(s-s_{\alpha1})(s-s_{\alpha2})}
$$

I know that this function does not have any zeros, because of the 1.

From the solutions I can see that the poles are $s_{\alpha1}=-1+j$ and $s_{\alpha2}=-1-j$ and to tell you the truth, I have no idea how to get solution for the poles.

I tried to solve it using quadratic equation, but the only thing I get is a big mess.

Thanks!!


texts: Can I apply Fourier Transform to a non-time-indexed signal?

texts:
Say I have a signal that is not x-indexed. That is, the x-axis of the signal is the distance traversed by car and the y-axis is the heading direction of the car at the corresponding distance.

Can I apply the Fourier Transform to this signal?
If so, what is the physical meaning of this transformation? I believe that the horizontal axis is no longer frequency any more. What is it in this case?