texts: How to calculate a significant value of a definite curve/shaded region in a graph?

texts:
The attached graph below consists of $4.5$ seconds of pattern which I have recorded via a hardware device.

The $x$-axis is the timestamp value in seconds and the $y$-axis is the hardware value. The values are taken on $0.01$ seconds of interval.

As seen in the graph, there is a definite pattern/curve from $0.5$ to $3$ seconds. I tried calculating a definite value using the below steps:

Taking 2.5 seconds block starting from $0$ Blocks: $0 - 2.5$, $0.01 - 2.51$, $0.02 - 2.52$, etc.
Now divide this $2.5$ second block into $5$ parts of $0.5$ seconds interval blocks. $(2.5/5 = 0.5)$
So, first block of $0 - 2.5$ will give: $0, 0.5, 1.0, 1.5, 2.0, 2.5$
Now I have calculated the slope of these values $(x_2-x_1/y_2-y_1)$. This will give me $5$ slope values.
Lastly I have added these values, $(\textrm{Slope}_1 + \textrm{Slope}_2+\ldots + \textrm{Slope}_5)/5$ to get a single average value.

The thing is this value is definite for the time block of $0.5 - 3$ seconds (the curve which I want to identify), but sometimes it comes for other time blocks also.

Also, second issue is, In the current graph, most of the values start from near $0$ or at $0$ on the $y$-axis. In real time, values can start from $-1$ or $1$ also, but the curve will still be the same, I have tested it.

I think the only way is to calculate the value of the curve, as the $y$ values are not definite, but the curve shape is. In real time I would have around $30$ seconds of data and I need to check if any $2.5$ seconds from those $30$ seconds contains a curve like this.


texts: Help with deconvolving photon reads distributed across neighboring pixels

texts:
I'm building a radiation detector that collects photons in CCD pixels and we can relate the energy of the photon from the intensity of the pixel. To test the detector, I took the following spectrum of Iron 55, which is expected to have a single peak at 6 keV:

The tail at lower values of the energy is interesting because it's not supposed to be there. Taking a closer look at pixels in this range, I see that they are clustered next to each other, with what seems to be energies from a single photon distributed into multiple neighboring pixels.

I am intrigued by the possibility of getting a cleaner spectrum by collecting these pixels together. I could try manually going through and clustering at each pixel, but I thought a more elegant way would be to deconvolve the image using a simple point spread function. In the result of this deconvolution I need 2 things to be true:

the energy of the system should be conserved, so the sum of the pixels (the $L_1$ norm) should be conserved.
there should be no negative pixel intensities, since those are unphysical.

I did that, using a Gaussian PSF with a half-pixel standard deviation. However, I find that I either need to violate condition 1 of my requirements, or condition 2 of my requirements. Just using the gaussian PSF with $L_1$ norm 1, which conserves the energy, gives me a plot like this:

The white pixels are negative, unphysical values. At this point I'm stumped, but also obsessed with this problem, because it seems like it should have a clean solution. Any advice?


texts:
I had a few questions on sampling(I'm quite new witht his), I tried to answer them, I think that I did the first one correct , but not sure about the 2 other:
.
given the next functions,Is it possible to discretely sample the function? if so what is the maximal allowed distance between the samples? otherwise explain why?

$$f(x)=\sin(\alpha x)$$

For this one I said it is possible using Nyquist theorem, assumin $\alpha =2\pi f$ and $T=\frac{1}{f}=\frac{2\pi}{\omega} $ then the allowed distance is $\frac{T}{2}=\frac{\pi}{\omega}$

(hopefully i got this allright).

now from the second question I'm not sure.

$$f(x) =
\left \{
\begin{array}{cl}
1,  & \text{$-1\leq x \leq 1$} \\
0, & \text{else}
\end{array}
\right .
$$

for this function I'm pretty sure the answer is it is not possible, since it is a straight line, so I'm not sure.. (please help me with this one).

$$f(x)= \mbox{convolution between function of question 1 , function of question 2.}$$

In this case I'm guessing it is possible since its the same as sampling first question function (only that it is sliced).


texts: how to detect number of syllables in a batch of audio recordings?

texts:
I have a corpus of short speech samples from Kiswahili speakers, and I want to detect the number of syllables in each audio recording. How should I approach this task?

Background: I asked a conceptual/programming variant of this question on StackOverflow showing my attempt at peak detection in R. Users migrated the question to CrossValidated, but someone suggested that SP might be a better home. I asked how to handle the situation on SP Meta, and a user gave me helpful advice to frame this question more conceptually on SP.


texts: Digital receivers for the reception of non pulse shaped signals

texts:
We are used to deal with pulse shaped signals, for example signals that go through root raised cosine filters at the transmit, and then  matched with root raised cosine filters at the receive

In addition to these functions there are carrier phase and frequency synchronization, symbol timing synchronization at the output of the matched filter , equalization,  etc..

Very briefly please, how would things be different at the receiver, if any, when receiving a signal that has not been pulse shaped by any conventional pulse shaping filter at the transmit(for whatever reason)?


texts:
I am recently acquiring data through ADC on signals ranging from $0.1\textrm{ Hz}$ to $50\textrm{ kHz}$ ($F_s = 128\textrm{ kHz}$). I noticed that apart from the white noise, I have something similar to $1/f$ noise with a cutoff frequency near $1\textrm{ Hz}$. There isn't much that I found on equalizing this $1/f$ noise on google. I am wondering what type of filters or algorithms that I can use to equalize the noise and have a flat noise bed in my frequency range as post-processing, so that I can see my input signals without the decay of magnitude caused by $1/f$ noise.

Also, I am wondering if I am able to equalize this noise, will it cause distortion in my signal in the time domain, or phase response?

Any advice, links to references or help would be greatly appreciated.


texts: taylor expansion of scale space function

texts:
I see the following expression from http://en.wikipedia.org/wiki/Scale-invariant_feature_transform

The quadratic Taylor expansion of the Difference-of-Gaussian scale-space function,  with the candidate keypoint as the origin is

$$D(\textbf{x}) = D + \frac{\partial D^T}{\partial \textbf{x}}\textbf{x} + \frac{1}{2}\textbf{x}^T \frac{\partial^2 D}{\partial \textbf{x}^2} \textbf{x}$$

where D and its derivatives are evaluated at the candidate keypoint and $\textbf x = (x,y,\sigma)$  is the offset from this point.

I'm very confused about this expression. I don't know why $D^T$ appears and don't understand everything. Is there someone to help me?


texts: What is zero-phase filtering and forward-backward filtering?

texts:
I am pretty familiar with IIR and FIR filtering and I have implemented them in my several projects. However, recently i found something called zero phase filtering. I have tried to understand it but not going well. All i know is that zero phase filtering is filtering that do not produce any phase delay and distortion

Is it right?

I have understood that FIR have a constant delay for a group signal while IIR have a non linear delay based on the frequency, so it will distort the signal phase.so there is a zero phase IIR to overcome the cons

How does the zero phase IIR work? (Needs an overview and maybe some examples)

And i have read about forward backward filtering that also remove the delay and distortion of filtering

How does it work?
Is it one of the zero phase filtering method?
It said that it can be done in offline filtering. What is offline filtering?

Thank you


texts: FFT doesn't always have frequency in all frames

texts:
I'm trying to do an FFT on live audio to see its spectrum.

So I have a sampling rate of 44100 Hz and from this take out 1024 bins on my FFT

So to test it I play a continuous tone at 15000 Hz and another one at 16000 Hz and look for it on my spectrum. I try to use higher ones to cut out noise from around. I have made the sound in a way there is no clipping

The thing is when I look at the individual frames of the FFT in each 2048 byte block of data I can see the peaks at these frequencies but some frames won't have the peaks, yet the others would.

How can this happen? How can I avoid it or what can I do to the way I'm reading the audio so that I can see the two tones in every frame?


texts:
I have a system that gives me a noisy data set similar to the one generated by this matlab/octave code. The y-axis represents the signal intensity and the x-axis represents spatial distance.

noise = randn(1,512);
points = gaussian(512,.01);
noisy = points + .02*(noise');
plot(noisy)

The actual signal will be asymmetric and not the result of a gaussian function(but close). I can't create the noise accurately by hand and the underlying signal algorithmically so image these two together.

The goal is to construct an interpolating curve, using bicubic splines, based on only one data set of this signal. The signal cannot be captured multiple times and averaged together to remove the noise. Since the signal is noisy, I can't simply pick points every N samples as the interpolation points so I need to do some signal conditioning first.

I will not be able to use an exact model of the underlying physical process but the sources of noise are known and I can make some assumptions regarding the maximum slope and SNR. The interpolation is used to correct for manufacturing limitations so the curve is intended to be smooth(exclusive of noise) but for cost reasons is not. Therefore we don't expect any large deviations from this example.

My initial approaches used an M-point moving average filter to smooth the signal, then used polyfit handle the ends of the signal. This worked but could potentially reduce the center peak, if M is large enough. I've also since learned smoothing signals before interpolation is generally discouraged.

What would be a better approach to achieve this?


texts: Ripple in time-domain using lanczos resampling

texts:
I'm seeing ripple in the time domain when resampling a "unity signal", that is, all ones, with lanczos filters. The ripple/error exactly resembles a cosine with a period of the sampling rate - ie. the error is largest at a fractionate sample offset of 0.5. Using larger $a$ for the lanczos kernel reduces the error, but it never ceases to exist.

Is this supposed to happen, or am I making an implementation error somewhere? Doing the convolution by hand shows that the terms never sum up to 1 (unless at offset = 0). Truncated sinc shows the same problem.
Are there any windowed sinc filters that do not exhibit this kind of error? Otherwise, any other suggestions for resampling techniques? I'm generally looking for something that offers the performance of a lanczos kernel with $a$ > 10.

Edit: Thanks for the comments. I'll elaborate a little on what I want: So I'm interpolating a waveform display. What's important and interesting is, that the interpolation passes through the original sample points and resembles the 'true' shape of the signal as much as possible. This means this is good:

So under- and overshoot is not a problem in itself (technically, it's the whole point). This, however, is (note this is just a DC signal) - ie. nothing in the signal provokes the ringing, it is just an artifact depending on your fractionate interpolation offset:


texts: Bandpass filter to get EEG frequency bands?

texts:
I'm doing a project with an EEG to detect driver drowsiness and learned about EEG frequency bands (alpha, beta, gamma...).

As the frequency bands are simple frequency ranges, I wonder if I can use several bandpass filters to get them (instead of using WPT / FFT)?

Is this approach correct? Is there any reason not to do it (performance)?


texts: Why unit impulse function is used to find impulse response of an LTI system?

texts: Hello i am working in digital image restoration field, recently i have studied concept of convolution, i studied that to find the impulse response/point-spread function of an LTI system, an unit impulse function is passed from the system and its output is considered as impulse response, My question is that, Why an unit impulse function is used to find the impulse response of any LTI system, what is the mathematical and practical reason behind it.


texts: Fastest Algorithm for Distance Transform

texts:
I am looking for the fastest available algorithm for distance transform.

According to this site http://homepages.inf.ed.ac.uk/rbf/HIPR2/distance.htm, it describes:

  The distance transform can be calculated much more efficiently using
  clever algorithms in only two passes (e.g. Rosenfeld and Pfaltz
  1968).

Searching around, I found: "Rosenfeld, A and Pfaltz, J L. 1968. Distance Functions on Digital Pictures. Pattern Recognition, 1, 33-61."

But I believe we should have a better and faster algorithm than the one in 1968 already? In fact, I could not find the source from 1968, so any help is highly appreciated.


texts: Distribution of a signal covariance matrix

texts:
A common estimation problem in signal processing assumes the following signal model
\begin{equation}
 \mathbf{r} = \sum_{i=1}^{Q}\alpha_i\mathbf{s}\left(w_i\right)+\mathbf{n}
\end{equation}
where $\mathbf{s}(w)$ is a known vector function, and $\mathbf{n}\sim \mathcal{CN}(0,\sigma^2\mathbf{I})$ is white Gaussian noise with variance $\sigma^2$.

Column vector $\mathbf{r}\in\mathbb{C}^N$ constitutes the noisy data.
Parameters $w_1,\ldots,w_Q$ are unknown and need to be estimated.
$\alpha_1,\ldots,\alpha_Q$ are zero-mean complex random Gaussian variables with unknown variance.
$Q$ is the number of parameters to be estimated and is known.

A very useful family of estimators rely on some properties of the correlation matrix which is
\begin{equation}
 \operatorname{E}\mathbf{r}\mathbf{r}^H =
 \mathbf{S}\left(\mathbf{w}\right)
 \mathbf{A}
 \mathbf{S}^H\left(\mathbf{w}\right)
 + \sigma^2\mathbf{I}
\end{equation}
where $\mathbf{A}=\operatorname{E}\,[\alpha_1,\ldots,\alpha_Q]^T[\alpha_1,\ldots,\alpha_Q]^*$, $\mathbf{S}(\mathbf{w})=[\mathbf{s}(w_1),\ldots,\mathbf{s}(w_Q)]$, $(\cdot)^T$ is the transpose operator, $(\cdot)^*$ is the conjugate operator, and $(\cdot)^H$ is the transpose and conjugate operator. In practice, the correlation matrix is estimated using the sample covariance matrix defined as
\begin{equation}
 \hat{\mathbf{R}}=\frac{1}{K}\sum_{k=1}^{K}\mathbf{r}_k\mathbf{r}_k^H
\end{equation}
where $\mathbf{r}_k$ are independent realizations of $\mathbf{r}$.

Denote $\lambda_1,\ldots,\lambda_N$ and $\mathbf{v}_1,\ldots,\mathbf{v}_N$ the eigenvalues and eigenvectors of $\hat{\mathbf{R}}$, respectively. The eigenvalues are ordered from largest to smallest (all are non-negative because $\hat{\mathbf{R}}$ can be shown to be a Gramian matrix). Then in such case it turns out that
\begin{equation}
 \begin{bmatrix}
  \mathbf{v}_1 \cdots \mathbf{v}_Q
 \end{bmatrix}
 \operatorname{diag}\left(\lambda_1-\sigma^2,\ldots,\lambda_Q-\sigma^2\right)
 \begin{bmatrix}
  \mathbf{v}_1 \cdots \mathbf{v}_Q
 \end{bmatrix}^H
 \rightarrow
 \mathbf{S}\left(\mathbf{w}\right)
 \mathbf{A}
 \mathbf{S}^H\left(\mathbf{w}\right)
\end{equation}
in probability when $K\rightarrow +\infty$.

Without going into further unnecessary detail my question is

Is it known the exact/asymptotic/approximated distribution of $\begin{bmatrix}
  \mathbf{v}_1 \cdots \mathbf{v}_Q
 \end{bmatrix}
 \operatorname{diag}\left(\lambda_1,\ldots,\lambda_Q\right)
 \begin{bmatrix}
  \mathbf{v}_1 \cdots \mathbf{v}_Q
 \end{bmatrix}^H$?

As far as I know the sample covariance matrix has a Wishart distribution, but I do not know what is the distribution of such signal covariance matrix.


texts: Efficient correlation of a low duty cycle training sequence

texts:
Is there a way to efficiently correlate a training sequence that is N samples long, framed at M samples where M >> N, with L occurrences of such frames (see below). For the pedants, the training sequence may or may not vary in the L subsequent frames, thus the longer correlation template.

The samples in between the subsequent training sequences would be unknown so they would be 0's in the correlation template. There is a lot of wasted processing with the multiplications/additions of all the 0's. I am wondering what optimizations can be made by taking advantage of such a signal structure.

Example:

       <- N ->
        _____
       |     |         ... repeated L-1 times ...
 ______|     |_______
 <------- M --------->


texts:
In digital communications systems we say that the received signal can be expressed as

$$y(t)= x(t)+ n(t)\tag1$$

if x(t) is the transmitted signal and $n(t)$ is the noise.

I am wondering whether a better model would be

$$y(t)= x(t-t_0)+ n(t)\tag2$$

To account for delay caused by channel. Is my logic correct?

If so then, what are the assumptions made when using (1) and what are we accounting for if we use (2)? ie what type of errors are we accounting for (delay spread, or time acquistion error...)

I am looking forward to hearing your responses and discussion on this. thanks


texts: Finding pitch from a wideband spectrogram

texts: We can use the narrowband spectrogram for finding the pitch of an utterance. But is it possible to find the pitch using the wideband spectrogram? If so, how would that be done?


texts:
This post follows this previous resolved post where I was trying to find the inverse Z-transform of a more simple filter output (that I used as an example to get the methodology). The present filter implements an additional integrator and I tried to use the same method. The problem is, I am not sure about how to deal with repeated poles.

$H(z)$ and $X(z)$ are given by:

$ H(z)= \frac{(1-z^{-K})(1-z^{-L})}{{(1-z^{-1})}^{2}}$
$ X(z) =\mathcal{Z}\{x[nT]\}= \mathcal{Z}\left\{Ae^{-\frac{nT}{\tau}}u[n]\right\} = A\dfrac{z}{z-e^{-\frac{T}{\tau}}}= A\dfrac{1}{1-e^{-\frac{T}{\tau}}z^{-1}} $

I want to express $y[n]$ as a function of the pulse amplitude $A$, $L$, $K$, $n$ (sample number), $T$(sampling period), and the time constant $\tau$; my goal being to retrieve the pulse amplitude $A$.

I used polynomial long division to develop $H(z)$ and i came up with something like this, assuming $K>L$:

$(1)\space H(z)=1+2z^{-1}+3z^{-2}+...+(L+1)z^{-L+1}+Lz^{-L}+Lz^{-L-1}+..+Lz^{-K+1}+(L-1)z^{-K}+(L-2)z^{-K-1}+...+z^{-K-L+2}$

Now I describe $Y(z)$ using partial fraction expansion: $Y(z) = H(z).X(z) = A.\left[\frac{(1-z^{-K})(1-z^{-L})}{{(1-z^{-1})}^{2}(1-e^{-\frac{T}{\tau}}z^{-1})}\right] = A.\left[\frac{B_0}{(1-e^{-\frac{T}{\tau}}z^{-1})} + \frac{B_1}{(1-z^{-1})} + \frac{B_2}{{(1-z^{-1})}^2}\right] $
$$\begin{align*}
& (1-z^{-K})(1-z^{-L})=B_0.{(1-z^{-1})}^{2}+ B_1.(1-z^{-1}).(1-e^{-\frac{T}{\tau}}z^{-1})+B_2.(1-e^{-\frac{T}{\tau}}z^{-1}) \\
& H(z)(1-z^{-1})^{2}= B_0.{(1-z^{-1})}^{2}+ B_1.(1-z^{-1}).(1-e^{-\frac{T}{\tau}}z^{-1})+B_2.(1-e^{-\frac{T}{\tau}}z^{-1})\\
\end{align*}$$
We obtain:
$$\begin{align*}
&B_0 = H(z) = 1+2z^{-1}+3z^{-2}+...+(L+1)z^{-L+1}+Lz^{-L}+Lz^{-L-1}+..+Lz^{-K+1}+(L-1)z^{-K}+(L-2)z^{-K-1}+...+z^{-K-L+2}\\
& B_2 = 0 \\
\end{align*}$$
Not sure about $B_1$:
$$\begin{align*}
&(2)\space B_1 = \frac{d}{dz^{-1}}\left[H(z).(1-z^{-1})^{2}\right]_{z^{-1}=1} = \frac{d}{dz^{-1}}\left[\frac{1-z^{-K}-z^{-L}+z^{-K-L}}{1-e^{-\frac{T}{\tau}}z^{-1}}\right]_{z^{-1}=1} = \left[\frac{(Kz^{-K}+Lz^{-L}-(L+K)z^{-K-L})(1-e^{-\frac{T}{\tau}}z^{-1})-(1-z^{-K}-z^{-L}+z^{-K-L})(e^{-\frac{T}{\tau}}z^{-2})}{(1-e^{-\frac{T}{\tau}}z^{-1})^{2}}\right]_{z^{-1}=1} \\
\end{align*}$$
If $B_0$ is correct (I am not sure), the inverse z-transform of $\frac{B_0}{(1-e^{-\frac{T}{\tau}}z^{-1})}$ can easily be retrieved. As for  $B_1$, I get:
$$\begin{align*}
&(3)\space \frac{B_1}{1-z^{-1}}=-\frac{(1-z^{-K})(1-z^{-L})e^{-\frac{T}{\tau}}z^{-2}}{(1-e^{-\frac{T}{\tau}}z^{-1})^{2}(1-z^{-1})}=-\frac{(1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1})e^{-\frac{T}{\tau}}z^{-2}}{(1-e^{-\frac{T}{\tau}}z^{-1})^{2}}{}
\end{align*}$$

Now according to this table http://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html, my expression looks similar to the follow one.

However I tried to apply it on a concrete example (using Excel), and it looks like I have missed something or got something wrong.

Any help would be greatly appreciated


texts:
I'm trying to extract / fill the walls in a floor plan, the walls are filled with a (cross-)hatch pattern at a more or less fixed angle.

I've managed to remove most of the "noise" by masking the DFT so I have both the +45 deg lines and -45deg lines in separate images.

The boundaries of the walls are quite obvious to the human observer but somehow I'm stuck at segementing those lines into groups, i.e connected components.

What I could think of is to apply a run-length encoding and then looking for subvectors of length 4 which fullfill a certain ratio of black and white. Or a cross-correlation of a pre-defined hatching pattern. The issue I have with these methods is that I'd have to introduce a parameter which has to be tuned for every floor plan, this is what I want to avoid.

So I'm asking you guys whether you could give me some input how to approach this problem.


texts: In papers or textbooks, I have found several different definitions of PSD. I think I understand the meaning of PSD and all these papers or textbooks agrees that PSD and Fourier Series share a proportional relationship. But still there are different formulas from P(f)=1/T(〖F(f)〗^2), P(f)=〖∆t〗^2/T(〖F(f)〗^2), P(f)=1/4∆f(〖F(f)〗^2), P(f)=(〖F(f)〗^2) or others. From what I learn, I think on the point of conception of power(or energy) and Parseval's theorem, discrete PSD should be described as simply P(f)=1/T(〖F(f)〗^2). So why is there so many different formulas? Are they wrong or do I misunderstand something?


texts:
I am generating a simple musical instrument using a ST board. basically I have a sensor which detects motion or spatial angle and generates sounds of differnt pitch and volume depending on the angle. I already have the driver for the audio codec so all I need to do is to generate sound samples and feed the samples to it. Now I am able to get the angle readings from the sensor but the challenging part is how to generate sound. From google search so far, I think I need to generate a sin function with frequency and amplitude(volume) set according to the reading from the sensor. But my codec assumes a sampling rate 0f 48KHz so how would I go about generating sinusoids with different frequencies for a fixed sampling rate?

So far I have done this:

samplingRate = 48000;
n = 0;//reset once there is a change in frequency

//this function is called 48000 times a second
int generateSineWave(float frequency,float volume)
{
  int temp = volume*(sin(2*pi*frequency*n) ;
  n = n + 1;
  if(n == samplingRate )
    n = 0;
  if(abs(temp) > MAXVAL)
    return ERROR_CODE;
  return temp;
}

This seems to be working (I am hearing something),but I am not sure if it's generating the right frequency sinusoid. also the sound I am hearing is not very pleasant, how would I go about generating complex tones(like the ones in piano for example)? I guess also my control variables(spatial angles) needs to be low pass filtered. But apart from that any idea on how I could generate more audibly pleasant waves?

I am a newbie in sound synthesis so any suggestion/help would be greatly appreciated.


texts: Image datasets used for salient objects?

texts: I would like to understand reasons for different image datasets used for salient object detection: their difficulties, their origins and likely future datasets. Does there exist any work explaining them such their progression and characteristics over time?


texts:
If we filter out ideal white noise using an ideal LPF of cutoff frequency 10 KHz and then sample it at 30 KHz , is the resulting discrete signal statistically independent? I would like to know the statistical behaviour of the output signal.

I was attending an on-line test. My answer was told wrong.Below is the snapshot from the on-line test.

I would also like to know what would be the result if sampling frequency is below Nyquist frequency?


texts: Is DCT (Discrete Cosine Transform) of Type-2 lossless or lossy?

texts:
As I know, in general, the DCT is lossless. But I'm not exactly sure about Type-2 of DCT.
Is it lossless or lossy?


texts: Cannot reconstruct constant overlap add STFT perfectly

texts:
I am implementing an STFT function for spectrogram of certain audio files.
Here's my function:

def stft( input_sound, dft_size, hop_size, zero_pad, window):

    if hop_size != dft_size:
        number_frames = int(np.ceil((len(input_sound) - hop_size) / (dft_size
        - hop_size)))
    else:
        number_frames = len(input_sound)//dft_size

    matrix = np.zeros(shape = (dft_size//2 +1, number_frames), dtype =
    complex)

    for i in range(number_frames):
        if hop_size != dft_size:
            block = input_sound[i*(dft_size-hop_size): i*(dft_size-
            hop_size)+dft_size]
        else:
            block = input_sound[i*(dft_size): i*(dft_size)+dft_size]

        #  pad zeros to have enough samples to compose the last frame
        if len(block) < dft_size:
            block = np.append(block, np.zeros(dft_size-len(block)))
        # apply window
        block = block*window
        # take dft
        dft = (np.fft.rfft(block, n=dft_size+zero_pad)) # dft_size/2 + 1

    matrix.T[i] = dft

    return matrix

Here's my inverse STFT function:

def istft( stft_output, dft_size, hop_size, zero_pad, window):

    if dft_size != hop_size:
        length = dft_size*len(stft_output.T) - hop_size*len(stft_output.T) +
         hop_size
    else:
        length = dft_size*len(stft_output.T)
    output_sound = np.zeros(length)

    for i in range(len(stft_output.T)):
        dft = stft_output.T[i]
        idft = np.fft.irfft(dft, n=dft_size+zero_pad)
        idft = idft*window
        if dft_size != hop_size:
            output_sound[i*(dft_size-hop_size): i*(dft_size-
            hop_size)+dft_size] += idft
        else:
            output_sound[i*(dft_size): i*(dft_size)+dft_size] += idft

    return output_sound

My Inverse STFT code works as they did reconstruct the sound back. However, my problem is that I can not get it reconstruct perfectly.
The audio file I am using right now is a drum clip and its sampling frequency is 11025 Hz with total 28884 samples. The parameters I used is dft_size = 1024, hop_size = 512, window = np.hanning(1024), zero_pad= none.

As I researched online, for hanning window, hop size being 1/2 of its windows length should get you perfectly reconstruction. But I didn't. So I wonder if there is something I coded wrong or I haven't found the right combination of dft size, hop size, and windows yet.


texts: Determine the Signal Curve from Parameters of a Power Curve by Noisy Measurement

texts:
I have a class of signals described by function:
$$
f(inc,d,t)=inc\cdot t^d
$$
where inc and d have a finite set of values like 1, 2, 3, i.e.
$$
inc, d\in \left \{1,2,3  \right \}
$$
and
$$
0\leq t<1
$$
Example plots:

I need to determine parameters inc and d given a discrete signal of f obscured by additive white Gaussian noise. Preferably an algorithm capable of working at real time.
Some candidates are:

sole differentiation (this amplifies noise)
filtering + differentiation (loses information)
curve fitting (least squares?)

Polynomial fitting seems to be an overkill, because function is known and is monotonic. Isotonic regression also seems too general, because it is applicable to all increasing/decreasing functions. Maybe there is a better solution when we have a concrete function and a set of possible parameters.


texts:
There is this following question.
Consider the transformation H {x}[n] = n x [n]. Does H define an LTI system?
What I understood from the question is:
x[n] -> H -> y[n]
or according to the question,
x[n] -> H -> n x[n]
to prove the time invariance, if I put n-no instead of n, it should be,
x[n-no] -> H -> y[n-no] = n x[n-no]
but it turns out to be,
x[n-no] -> H -> n x[n - no] - no x[n-no]
which means the output depends on extra term no x[n-n0] and is not time invariant.
Incase, the trasformation was any constant * x[n] the it would have conserve the time invariance property, if I am not wrong.
But the proof mentioned in the notes is,

I have difficulty understanding the proof. If anyone could explain it to me or if I understood the question wrong, I'd appreciate it.


texts: Detection of a Circle in Noisy Image Data

texts:
I have an image that looks like the one below:

I'm trying to find the radius (or diameter) of the circle.
I have tried using circular Hough transform (via matlab's imfindcircles(bw,[rmin rmax],'ObjectPolarity','bright')) , and by fitting to a circle or an ellipse (home made function that works pretty well for less noisy data, see below).

I've also tried some image processing to get a clearer circle, for example, see below:

se = strel('disk', 2);
bw = imdilate(bw, se);
bw = bwareaopen(bw,100000);
bw =  edge(bw);

However, when I feed the processed image to either techniques (Hough and circle\ellipse fitting) neither of them manage to detect the circle in a decent manner.

Here's a code snippet of the circle finder I wrote (matlab)
       [row col]=find(bw);
       contour = bwtraceboundary(bw, row(1), col(1)], 'N', connectivity, num_points);

    x = contour(:,2);
    y = contour(:,1);

    % solve for parameters a, b, and c in the least-squares sense by
    % using the backslash operator
    abc = [x y ones(length(x),1)] \ -(x.^2+y.^2);
    a = abc(1); b = abc(2); c = abc(3);

    % calculate the location of the center and the radius
    xc = -a/2;
    yc = -b/2;
    radius  =  sqrt((xc^2+yc^2)-c);

Alternative approaches will be appreciated...


texts: How should I modify my filter coefficients for the filter to be unaltered, when I change my sample rate

texts:
I am guessing this is a commonly encountered problem, but I am lost on how to really solve it.

I have FIR filter coefficients that are used to match the frequency response of data sampled at 3 kHz. I have 12 coefficients, which I would like to implement as a filter for data being sampled at 48 kHz, but I understand, I can't directly do this.

So far, I've tried this (MATLAB):

Find roots of the 12 order polynomial (which are the zeros).
Find the magnitude and angle of the zeros.
Scale the magnitude(m) and angle(phi) by a factor of 3 kHz / 48 kHz.
Construct new zeros by multiplying m .* exp(phi);
Find the polynomial whose roots are the new zeros using poly.
Coefficients of the new polynomial are the new filter?

However, I am not sure how sound this logic really is, or if there's something wrong here entirely. I also think I'm losing a lot of information to resolution.

Any pointers on how to go about this would be really helpful!


texts: The spectrum of the sum of two periodic signals

texts:
Considering the continue, periodic signals $s_1(t)$ and $s_2(t)$, with the period $P_1$ and $P_2$ respectively.

Consider now a new signal $s$ the sum of the two pervious signals.

I come to know that the spectrum of the sum, is not necessarily the sum of the spectrums unless the two signals $s_1(t)$ and $s_2(t)$ have the same period. With spectrum I mean:

$$X\in [\mathbb{Z} \rightarrow \mathbb{C} ]$$

$$X[k] = \langle x(t), \phi_k\rangle = \frac{1}{P} \int_P x(t)e^{-ik\omega_0 t}$$

with $\omega_0 = \frac{2\pi}{P}$. $P$ here represents the period.

Note that $x(t)\in [\mathbb{R} \rightarrow \mathbb{C}]$ is periodic with the period $P$.

I don't understand why? Because at it seems this violates the linearity properties.


texts:
I am programming the simplest way of implementing and IIR filter. However, I am struggling with the use of buffers. Do any of you have any simple implementation of an IIR filter implementing manually the convolution? I would like to have something simple, implemented on a function or something like that on C/C++.
thanks in advance.


texts: Spectral density interpolation: Applications?

texts:
Let $\Phi_1$ and $\Phi_2$ be two matrix-valued power spectral densities. I wonder whether the problem of interpolating $\Phi_1$ and $\Phi_2$ has been studied in the literature and/or applied in practical scenarios.

In case of 1-D spectra, I know that a possible application could be speech/voice morphing, however I cannot imagine a practical application in which the considered spectra are matrix-valued functions.

Any suggestion/pointer to the literature is really appreciated!

  Addendum. The motivation behind my question comes from the fact that I'm currently working on the geometry of the space of power spectral densities (mostly from a mathematical viewpoint) and I would like to apply some of the mathematical tools I analyzed to some "practical" situations. One of these tools is a method for interpolating spectral densities. I'm especially interested in the case in which spectral densities are functions that take matrix values, since in this case the geometry of the space of spectral densities becomes particularly "rich" and "interesting".


texts: Mirrored spectrum after computing derivatives of 2D function in Fourier domain

texts:
I use FFTW library to compute gradients and divergences of my function by the relation
$\nabla_\vec x = i \vec k $. I use the transformation functions for the real signal to save the memory and indeed my signal is real. Initial conditions are set so that only low frequency disturbance is present and for what I expect it should remain in the lower part of the spectrum. Nevertheless, after dozens of computations the spectrum becomes mirrored to high frequencies as on the attached picture. Looks like Van Gogh painting but it's not what I actually need:).
I suspect that this is some numerical artefact. If you've ever seen something like this please let me know what could that be. On the axis there are numbers of samples and the 0th frequencies are in the middle.


texts: Why is my filter unstable and self-oscillating in this case?

texts:
I am using resonant bandpasses to simulate modes of guitar/piano strings for audio synthesis. While attempting to introduce coupling effects I have encountered a problem.

I will illustrate with some simple code:

exciterFinal = exciter - (output * scalar_0to1);

resoBP.setInput(exciterFinal);
resoBP.setFilterCutoff(frequency);
resoBP.setFilterQ(Q);

output = resoBP.getOutput();

In other words, I am subtracting the prior sample's output multiplied by a scalar from 0 to 1 from the usual input of the bandpass to simulate a damping effect.

The most extreme effect would be with a scalar of 1 which should immediately damp the bandpass completely.

This works fine at frequencies below around ~9.4 kHz with normal damping behavior at 48 kHz sample rate. But once you go above that frequency with a scalar of 1, the filters tend to immediately self oscillate to overload. I tried with a State Variable Filter and a Butterworth Filter and had the same result with both.

Increasing the sample rate to 96 kHz allows it to work up to about 17.7 kHz. Above that point it still oscillates. So sampling seems to be the limiting factor. But I cannot afford the processing power to go up this high or higher to get the full spectrum of sound.

What is the reason for this behavior and why does it only happen at frequencies above roughly 19% of the sample rate? Is there any way to fix it?

Update: It just seems the filters become unstable when they are too aggressively damped this way at high frequency and low sample rate. My solution is to damp the filters instead by modifying their Q values when I need aggressive damping and be more gentle when I need to use the subtraction-from-input method. Good enough I guess.


texts: Convolving time reversed discrete input with time reversed discrete impulse response

texts:
Question 3.1 d from Chapter 3 of Oppenheim's Signals & Systems
I have

$$x[n]=(-1)^n{u[-n]-u[-n-8]}$$
$$h[n]=u[n]-u[n-8]$$

and the question asks for $y[n]$ which is the convolution of $x[n]$ and $h[n]$ namely

$$y[n]=x[n]\star h[n]$$

I'm stuck with the limits on the sum

$$\sum_{k=-\infty}^{\infty}{\left((-1)^{k}(u[-k]-u[-k-8])\right)\left(u[n-k]-u[(n-8)-k]\right)}$$

Should the limit be '$0$' or '$n-k$' or '$n-8$'?


texts: How to detect a spiky signal embedded in a periodic signal

texts:
Say I've got a signal that I know is composed of periodic components. However, I'm worried that there are "spiky" events that are embedded in this signal, and potentially occurring periodically. I can't just look at the power spectrum of the signal, because the spiky event will show up across many frequencies in the power spectrum.

What kinds of methods are out there to tease apart components of a signal that can reasonably be decomposed into sine waves, from those that cannot (like a spike)?


texts:
I'm trying to get Parseval theorem working on a one-sided FFT. So far I have this code (matlab):

%% ODD
a = [1 2 3 4 5];
A = fft(a);
A1 = abs(A);
A2 = abs([A(1) 2*A(2:3)]);

Ea=sum(a.^2)
EA1=sum(A1.^2)/5
EA2=sum(A2.^2)/5

%% EVEN
b  = [1 2 3 4 5 6];
B  = fft(b);
B1 = abs(B);
B2 = abs([B(1) 2*B(2:3) B(4)]);

Eb =sum(b.^2)
EB1=sum(B1.^2)/6
EB2=sum(B2.^2)/6

And the output is:

Ea =  55
EA1 =  55
EA2 =  65
Eb =  91
EB1 =  91
EB2 =  107

The scaling for the one-sided is to preserve amplitude information and not for calculating Parseval theorem?


texts: Audio Signal Separation - Identifying Interest points

texts:
I asked this question here: Audio Separation of .wav signal but it wasn't clear, so, here is my second attempt:

First off, assume that I have a .wav file containing a sentence as follows:

"My name is Michael" I would like to extract, from this, the following:

"My" -> Phoneme (1)

"Name" -> Phoneme (2)

"Is" -> Phoneme (3)

"Michael" -> Phoneme (4)

This means that I have taken my 1D signal, and split it into a 2D signal (vector) that contains these particular words/phonemes which I can then analyse and identify. I would therefore, like to compute this in the time domain and not the frequency domain. Just to clarify again:

I take in a 1D signal containing a sentence, split this sentence up into different parts which contain this data: vect[0], vect[1], .... vect[4] Let's say in matlab I did the following command wavwrite(vect[0], ....) then it would output the word "My" and putting all the blocks together would give me the full sentence back.

Here is my "real-world problem" instead of Phonemes, I have bat calls, the length of each bat call is unknown at this stage. But here is a typical sample of a bat call: Here and for each of these bat calls, these need to be separated from the inputted signal and stored inside a vector (Just like the example above), this, then allows me to identify each of the bats and perform analysis on them.

This, like, the sample would give me a 2D vector containing each of the bat calls: "Bat1", "Bat2" ..... "Bat[n]" it is unknown the amount of time the bats have been recorded for, or, what the length of each of the bat call therefore is.

What I have done so far:

I have obtained the bat signal, processed it and I am given the following (Which is plotted):

I have also Emphasised the signal using the following formula:

rawSignal[i] = rawSignal[i] - (0.95 * rawSignal[i-1]);

And then I have Compressed the signal using the following:

    float param = 1.0;
for(unsigned i=0; (i < rawSignal.size()); i++)
{
    int sign = getSignOf(rawSignal[i]);
    float norm = normalise_abs_value(rawSignal[i]);
    norm = 1.0 - pow(1.0 - norm, param);
    rawSignal[i] = denormalize_value(norm, sign);
}

Which then gives me an output of the following:

I'm unclear to where I should go from here in identifying single elements ("calls") from this signal. Since, if I use zero-crossing and/or calculating the total energy of the signal and thus using a threshold then it will just remove the noise and I'm left with a compressed version of the signal.

Speaking to someone, they suggested that I should try and use the Cochleagram domain however, I'm not familiar with this and there is very little research on this available.

If anyone has any suggestions, or the algorithms that I could use then please suggest them.


texts:
I come from Computer Science so please pardon for my possibly wrong terminology.

I need to design a filter which has coefficients
$$h_0, h_1, \ldots, h_n, \ldots \quad\text{such that}\quad h_0 > h_1> \ldots > h_n > \ldots$$

The input is $x_n$. The output is $y_n$. The relationship between $h_n$, $x_n$ and $y_n$ is a convolution:
$$y_n = h_0 x_n + h_1 x_{n-1} + \ldots + h_n x_0$$

The problem is that I have to carry out a convolution to compute the output $y_n$ upon the arrival of every new input $x_n$. So I have to keep all $x_1, x_2, \ldots, x_n, \ldots$ in memory and do $n+1$ multiplications and $n$ additions every time.

I would like to ask whether there is any way to design the coefficient $h_n$ such that I can compute $y_n$ more efficiently. I know of one solution, which is $h_n=\lambda^n$ with $0<\lambda<1$ but this function decays too fast.

Any idea to approximate this procedure is also welcomed.


texts: Convert Fourier transform of a signal to $\omega$ form?

texts:
Suppose we have this signal:

$$x(t)=5 \sin(2\pi 1000 t) \cos(2\pi 10000 t)$$

I calculated the Fourier transform as below:

$$X(f)=\left(j\frac{5}{4}\right)\left[\delta(f-9000)+\delta(f+11000)-\delta(f-11000)-\delta(f+9000)\right]$$

and then I would like to convert it to $\omega$ form. I know that the result is

$$X(\omega )=\left(j\frac{5}{2}\right)\left[\delta(\omega -9000)+\delta(\omega +11000)-\delta(\omega -11000)-\delta(\omega +9000)\right]$$

I know that
$\omega=2\pi f T$
where $T$ is sampling rate.

What I cannot understand is how we convert the Fourier transform of our signal to $\omega$ form using above formula.


texts: Does zero-forcing equalizer need known channel impulse response?

texts: I'm studying some of the basic equalizing structres and I understand how Zero-forcing works, but it seems to me that a known channel impulse response is needed. Am I right? If so, what's the point? I mean, you're not likely going to know how the channel is, so how Zero-forcing is in any way useful?


texts: Fractional Delay Filters and Cutoff Frequences

texts:
I am trying to implement an interpolator for arbitrary sampling rate conversion of a one-dimensional signal (a fractional delay filter). I am aware of the fact that interpolation is in general a problem of low-pass filtering, and impulse response and passband widths are inversely proportionate (for the same kernel shape). However, I'm a bit uncertain of how precisely to apply this to fractional delay filters for different interpolation rates.

As I understand it, the frequency cutoff for the interpolation kernel should be at pi*low_rate/high_rate radians. But for what class of interpolation filters is this the case? E.g. for cubic B-spline interpolation (a non-Nyquist(1) piecewise polynomial, where Nyquist(1) means the filter's impulse response is delta when the source and target sampling rates are identical) for supersampling the kernel is based purely on the fractional position, and does not have any kind of provision for a cutoff, despite it being possible to compute an impulse response for the interpolator (which would seem to make the kernel frequency-dependent by the impulse response/passband duality mentioned previously). What about for other classes of interpolators such as linear (Nyquist(1), polynomial), sinc (Nyquist(1), non-polynomial) or Gaussian (non-Nyquist(1), non-polynomial)?

Additionally, non-Nyquist(1) interpolators can by converted to Nyquist(1) with the use of a fixed prefilter, e.g. as described in Variable Sample Rate Conversion Techniques for the Advanced Receiver. If in fact a given interpolator must have its impulse response dilated/contracted based on cutoff frequency, this means there are now two stages for these non-Nyquist(1) interpolators: the fixed prefilter and the potentially variable interpolation. At what stage(s) would the correction for cutoff frequency be applied in such a system, and would it follow the same formula stated above? I would have expected frequency correction to be necessary in the prefilter, but the paper doesn't seem to indicate the use of frequency-dependent prefilters with cubic B-spline interpolation despite that interpolator, as already stated, having no accommodation for cutoff frequency in the core algorithm.


texts: MATLAB's $\tt dsp.ParametricEQFilter$ bandwidth question

texts:
Consider the following code:

Fs = 48000;
h = dsp.ParametricEQFilter('CenterFrequency', 1000, ...
                           'Bandwidth', 1400, ...
                           'PeakGaindB', 10, ...
                           'SampleRate', Fs);

fvtool(h,'Fs',Fs,'FrequencyScale','Log');

Let the two bandwidth $-3\rm dB$ points $f_1$ and $f_2$.
Therefore: $f_1 = 300\ \mathrm{Hz}, \ f_2 = 1700\ \textrm{Hz}$

However, in the graph I see the $-3\ \mathrm{dB}$ point (i.e. the $+7\ \rm dB$ gain) at $487\ \rm Hz$. What am I doing wrong?


texts: What is the conceptual difference between Kolmogorov-Sinai entropy, Shannon entropy and Boltzmann entropy? Are they interchangeable and mean the same? Where can I find a good lucid explanation about Kolmogorov entropy and why it is needed and its significance. Thank you.


texts:
I want to go from this image

into this one:

So basically I need to scale the white square. The authors of the paper claim that this can be done in four steps:

zero-padding in real space (image is treated as 1D vector).
discrete Fourier transform.
cropping the central part of the spectrum in frequency domain.
inverse Fourier transform.

I don't see how it is being done in practice and if this works. Can you help out?


texts: Finding an Image's Orientation With DFT Frequency Amplitudes

texts:
I'm looking to get some general idea of an image's orientation so that I can then rotate it to the nearest 90 degree angle.

My idea on how to do this is to take the DFT on the image, then do a "radar sweep" to find which angle has the greatest sum of frequency amplitudes along the line, and call the winning angle the orientation.

I scan from 0 to 180 degrees instead of from 0 to 360, and make the radar sweep double sided since negative and positive frequencies seem to belong together.

For example, if i take the DFT of this image:

The frequency amplitudes are:

And if I do a radar sweep, I'd find the winning line / angle here:

Using that information, I can know where the nearest 90 degrees is for that blue line, figure out how much i'd need to rotate the dft image to make that blue line vertical, but instead rotate the source image by that amount, to make that image orient towards the nearest 90 degrees.

Something like this, but actually level (I'm working out some bugs at the moment, in my implementation)

Is this method sound?  Are there any complications that I'm going to encounter when moving to more complex images such as complex patterns or photographs?

Edit: In the above I didn't explain the reasoning for what I am doing.

If you take the dft of simple images - like a single line, or a repeating pattern of black and white stripes - the resulting frequency amplitudes image has data that runs perpendicular to the direction of that line or those stripes.

Line this:

DFT'ing to this:

May (possibly naive!) understanding is that I should then be able to tell which direction (from the center) has the highest sum, and that will tell me the angle perpendicular to the largest represented orientation of the image.

From there, i can find how much i'd have to rotate that line to get to the nearest 90 degrees, and apply that rotation to the source image, to make it line up to it's nearest 90 degrees.

Is my thinking too simplistic, or is this a solid approach?


texts: How would addition with a constant in time domain affect my signal in frequency domain?

texts:
I have an audio signals from two microphones, and I have to make them equal for equal input.
I.e. calibrate the output of microphones to give the same resulting signal on both of them.
The problem is with amplitudes. Both sensors aren't equaly sensitive (logically. But I need the output to be.

The naive solution of this problem would be to find a constant difference between signal from mic A and mic B. And simply add it to the signal that requires that boost.

But, well, I am introducing deviation to signal. I am pushing it up along Y axces. This increases the top amplitude, but I am increasing the silence too.

Well, signal will not really change, on speaker will be heard the same, as silence line continues to be the line,

but I don't see what other implications I may have if I do this.

How would this affect the frequency spectrum for instance.

Note: I cannot do:

(x/max(x)*constant

normalization because I loose the info about initial amplitude of the signal. My measurement uses amplitudes and they must stay as they are, except that microphones should be calibrated.


texts:
I need to ask about insertion of Zero Symbols after QAM modulation in a OFDM system. Such that the 4 QAM has a constellation with 4 corner and 1 center symbol indicating zero. What advantages or disadvantages woild that offer ?

Is this done for guard bits or as pilot or what ?

So I need to add more information to the question. I just had a tutorial on OFDM system on MATLAB and what my Tutor told was that in addition to the Cyclic Prefix, we may also include these zero symbols after we have modulated the bits with QAM. Which appear as Zero Points in Constellation. He was saying that they have the purpose of controlling the ISI as well. But I was confused because ISI was controlled by Cyclic Prefix. ?

Also I am not sure about what effects this would have on OFDM Spectrum.


texts: How to evaluate performance of an ARMA, MA or AR model?

texts:
How to evaluate performance of a model after estimating ARMA/MA/AR parameters for any process x(n)?

How to regenerate back a process after estimating average parameters?
 what kind of performance measure can be used to evaluate the corresponding model?


texts:
I found a filter for Hankel Transform . As far as i understand this is effectively used as an approximation for $J_0(x)$ and $J_1(x)$.

But how will the change in input affect these filter terms (e.g. $J_0(r) \rightarrow J_0(2r)$)

References : http://onlinelibrary.wiley.com/doi/10.1046/j.1365-2478.1997.500292.x/full


texts: How to make frequency axis for even and odd FFT length?

texts: Can I get help on how to make frequency axis going from negative frequency to positive frequency, (in Hertz), that will be the x-axis in an FFT result, but given either an even length FFT, or odd length FFT. I am having some troubles making it in MATLAB. (Assume you know the sampling frequency f_s).


texts: How do I deal with occlusion, warping, or perspective distortion for object detection

texts:
I am using SIFT features for detecting logos in clothes. SIFT is quite powerful, as it is scale invariant, but how do I deal with the problem of occlusion and warping?

More specifically, the logos in clothes might be partially occluded, and can also be warped when the clothes are stretched or squeezed.

Any suggestion and directions? Thanks in advance!

Some example are showed here:

Logo on leg:


texts: Deriving finite impulse response for polylogarithm

texts:
As a part of my research i have to use the following z-transform in matlab 'filter' function so as to derive the convoluted signal from the original one.

$$\frac{1}{{\rm Li}_{k}(z^{-1}e^{-b})}$$

But since i read that the polylogarithm  can be expressed as a function only for specific values of k (k can take many values, not necessarily integers). Can anyone guide me as to how can i use it for my matlab codes. Will i have to derive the FIR. Can it be derived in this case ?


texts:
I'm studying image warping and I'm unsure how to warp a starting rectangular image into another image where the box corner points have been changed as I want:

(image from http://www.inf.ethz.ch/personal/pomarc/courses/CompPhoto/cpv03.pdf)

As far as I've understood, I need to:

1) Decide the four points (e.g. by using direct mapping on the border cases of the transformation)

2) Create a lattice of integer points on it (how?)

3) Do reverse mapping on those integer points to get the values in the original image

4) ??

I can't say I understood a lot of that, can somebody help me out please?


texts: Capture legitimate huge increase/drop in fuel level from noisy measurements

texts:
I have a fuel sensor installed in a lorry fuel tank. The sensor measures the current fuel level (in percentage) of the tank. However, due to the fact that the fuel is sloshing in the tank, the fuel level fluctuates wildly — sometimes more than 10%.

I decided to use Kalman filter to predict the actual fuel level, using the following model:

F(k+1) = F(k) - C
C = fuel consumed between k and k+1, which is estimated using
    the distance travelled between k and k+1 (measured using GPS),
    and the average fuel efficiency (fixed) of the lorry (e.g. 35 L/100km).

So far, Kalman filter works perfectly when the lorry is moving and consuming its fuel. However, it fails when 2 events happened that are not captured by the model:

refuel (large increase in fuel level in a short span of time)
fuel theft (large decrease in fuel level in a short span of time)

You can see an example below, on the left side of the graph, where a huge drop in fuel level (most probably due to fuel theft) and a huge increase in fuel level (due to refuel) are smoothed out by Kalman filter, thus the information are lost. On the right side of the graph, where the lorry is moving and consuming fuel, the prediction works well.

So my question is, how do I smooth the fluctuations of the fuel level, and at the same time, capture refuel and fuel theft events? Any ideas? Maybe use something else other than Kalman filter?


texts:
I am trying to implement a sobel filter to locate horizontal edges in a 160x120 pixel image that is captured from a OV7725 CMOS camera module that outputs RGB565 data.

I've done the RGB565 -> 8-bit grayscale conversion as can be seen in the first picture

However, when I try to do my sobel operation on the grayscale picture it gets a bit funky and the output is multiplied for some reason. I've been looking at my implementation but can't seem to find out what is wrong...

Below is my implementation in C and I don't really care about the pixel at the edges.

void sobelFilter(int8_t *inputImg, int8_t *outputImg)
{
    uint16_t row;
    uint16_t col;
    uint16_t pos;

    int16_t pixel_x;
    int16_t pixel_y;

    const int8_t sobel_h[3][3] = {
            {-1, 0, 1},
            {-2, 0, 2},
            {-1, 0, 1}
    };

    for (row = 1; row < 120 - 1; row++)
    {
        for (col = 1; col < 160 - 1; col++)
        {
            pos = 120 * row + col;

            pixel_x = inputImg[pos - 121] * sobel_h[0][0];
            pixel_x += inputImg[pos - 119] * sobel_h[0][2];

            pixel_x += inputImg[pos - 1] * sobel_h[1][0];
            pixel_x += inputImg[pos + 1] * sobel_h[1][2];

            pixel_x += inputImg[pos + 119] * sobel_h[2][0];
            pixel_x += inputImg[pos + 121] * sobel_h[2][2];

            outputImg[pos] = sqrt( (pixel_x * pixel_x) );) );
        }
    }
}

Is there something funky with my implementation that I've missed or am doing wrong? Like why is it producing 3 copies of the "object" and why is the filtered image not being filtered the whole way down?


texts: Simple analog anti aliasing filter good enough if signal frequency is far from sampling frequency ? What about equivalent time sampling?

texts:
I have one and only one frequency signal, a sine signal at f = 10 kHz that is sampled at Fs = 320 kHz.

I have read about applying a simple analog anti aliasing (AAA) filter then oversample, do some digital anti aliasing (DAA) and decimate to get the correct signal frequency.

I understand how this is useful if signals of interest are close to Fs/2 as it might be hard to find an AAA filter that has sharp enough roll-off, low cost, simplicity etc. And therefore DAA filter in combination with a simple AAA filter is often preferred.

But in my case I don't have a signal close to Fs/2 at all. So I was thinking maybe I can get a slower/cheaper ADC and do no oversampling, since I might be able to properly attenuate frequencies that might alias even using a simple AAA? Is what I am thinking correct?

Additional related question:
How would this work if I was to apply equivalent time sampling (ETS) technique? This means that the signal is 10 kHz and the sampling frequency appears to be 320 kHz as well, but 32 samples are taken over 32 periods of the 10 kHz signal to produce 1 period. So the sampling rate is actually 10 kHz + (10 kHz / 32) = 10.3125 kHz. What I am not able to wrap my head around: will Fs = 10.3125 kHz or 320 kHz? Will the same filter work for both cases? I assume Fs = 10.3125 kHz thus giving a problem as the transition area will be relatively very small: Fs - f = 3125 Hz.

Other thoughts:
I am using a DDS to generate the sine signal, it should be quite pure and maybe AAA is not "needed" / very important?


texts:
I am working with Walsh coefficients.

I know the intuitive understanding is almost that that they are the degree of connectivity, but it is there a better way of thinking about it? What is the definition of a Walsh coefficient besides "What you get when you take the Walsh transform".

Thanks


texts: Intertial navigation on android phone with Kalman filter

texts:
Ok, so suppose I've got a phone with gyroscope, compass and 3-d accelerometer.
I wanted to track position of the moving phone for about 1 minute with let's say 50 mm accuracy. Actually let's say that I wan't to create some sort of pedometer.
Is this possible?

I already know that the signal from accelerometer is too noisy to calculate position through double integration after really small amounts of time because of the giant drift.

But I want to know if it's possible with some domain knowledge using Kalman filter?
So in my case I know that the phone would be in the users hand and it wouldn't accelerate more than some amount X m/s^2 and that velocity wouldn't be more than some V m/s and so on. Would that be possible in that case?

Also how to create covariation matrix for the accelerometer, gyroscope and compass?


texts: Frequency response of a system given its block diagram

texts:
Given the following block diagram, I am asked to find the frequency response $H(f)$ of the system.

This is what I have done: The output of the first block is $j2\pi fX(f)$, and after going through $t_0$, we have $j2\pi fX(f)e^{-j2\pi t_0 f}$. Then, the input of the last block is $j2\pi fX(f)(1-e^{-j2\pi t_0 f})$ and the resulting output is $\frac{j2\pi fX(f)(1-e^{-j2\pi t_0 f})}{j2\pi f}=X(f)(1-e^{-j2\pi t_0 f})=Y(f)$. Finally, we have that $H(f)=\frac{Y(f)}{X(f)}=1-e^{-j2\pi t_0 f}$.

I know that my procedure can be written in a more formal way but, in general terms, is it correct? And if I was asked to plot $|H(f)|$, one way could be computing the square root of $|H(f)|^2=(1-e^{-j2\pi t_0 f})(1-e^{j2\pi t_0 f})=2-2\cos(2\pi t_0f)$, am I wrong?


texts:
I am trying to learn about signal processing and fft's applications within it. My understanding is that a classical use case for fft is to remove "noise" (ex. high frequencies that are not supposed to be in the original). However, while on StackOverflow, I came across this post that seems to suggest that this is a bad idea.

Am I misunderstanding the post or fft's use cases? If so, how would it be used?


texts: Smoothing for damped wave signal (fixed variance noise, but changing SNR)

texts:
My colleagues and I are arguing about how to smooth a damped wave signal.

This signal is corrupted by white noise of a steady magnitude.

However the signal damps out as it goes along. So, the SNR and the noise as a standard dev of the signal range, both decrease.

My colleagues say I can denoise it the same the entire time, because the magnitude of the noise is the same. I say I have to adapt, because the SNR is the guiding factor.

Who's right?


texts: What does an ADC DNL specification mean, exactly?

texts:
I'm looking at the differential non-linearity specification for an Analog-Digital converter. The spec sheet claims that the DNL ranges from -1LSB to +1LSB, with a typical value of $\pm$0.6LSB. What exactly does that mean? Is this some value that's taken statistically? Does this mean the 3$\sigma$ value is 0.6LSB? I.E., for any given code, the range of voltages over which it has jurisdiction is 99.7% likely to be no greater than 1.6LSB or less than 0.4LSB?

Suppose I wanted to model a 16-bit ADC with these specifications. How would I go about injecting DNL into the output codes? The spec sheet is completely unclear as to what 0.6LSB means. Is there some statistical convention? Or would it just be some random number which caps at $\pm$1LSB for every single code?


texts: Discrete Cosine Transform (DCT) Coefficient Distribution

texts:
I have two images :

Original Image

Binarize Image

I have applied Discrete Cosine Transform to the two images by dividing the 256x256 image into 8x8 blocks. After, I want to compare their DCT Coefficient Distributions.

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import numpy as np
import os.path
import scipy
import statistics

from numpy import pi
from numpy import sin
from numpy import zeros
from numpy import r_
from PIL import Image
from scipy.fftpack import fft, dct
from scipy import signal
from scipy import misc

if __name__ == '__main__':
    image_counter = 1

    # Opens the noisy image.
    noise_image_path = 'noise_images/' + str(image_counter) + '.png'
    noise_image = Image.open(noise_image_path)

    # Opens the binarize image
    ground_truth_image_path = 'ground_truth_noise_patches/' + str(image_counter) + '.png'
    ground_truth_image = Image.open( ground_truth_image_path)

    #Converts the images into Ndarray
    noise_image = np.array(noise_image)
    ground_truth_image = np.array(ground_truth_image)

    #Create variables `noise_dct_data` and `ground_truth_dct_data` where the DCT coefficients of the two images will be stored.
    noise_image_size = noise_image.shape
    noise_dct_data = np.zeros(noise_image_size)
    ground_truth_image_size = ground_truth_image.shape
    ground_truth_dct_data = np.zeros(ground_truth_image_size)

    for i in r_[:noise_image_size[0]:8]:
         for j in r_[:noise_image_size[1]:8]:
              # Apply DCT to the two images every 8x8 block of it.
              noise_dct_data[i:(i+8),j:(j+8)] = dct(noise_image[i:(i+8),j:(j+8)])
              # Apply DCT to the binarize image every 8x8 block of it.
              ground_truth_dct_data[i:(i+8),j:(j+8)] = dct(ground_truth_image[i:(i+8),j:(j+8)])

The above code gets the DCT of the two images. I want to create their DCT Coefficient Distribution just like the image below:

The thing is I dont know how to plot it. Below is what I did:

#Convert 2D array to 1D array
noise_dct_data = noise_dct_data.ravel()
ground_truth_dct_data = ground_truth_dct_data.ravel()

#I just used a Histogram!
n, bins, patches = plt.hist(ground_truth_dct_data, 2000, facecolor='blue', alpha=0.5)
plt.show()

n, bins, patches = plt.hist(noise_dct_data, 2000, facecolor='blue', alpha=0.5)
plt.show()

image_counter = image_counter + 1

My questions are :

What does the X and Y-axis in the figure represents?
Are the value stored in noise_dct_data and ground_truth_dct_data, the DCT coefficients?
Does the Y-axis represents the frequency of its corresponding DCT coefficients?
Is the histogram appropriate to represent the DCT coefficient distribution?
The DCT coefficients are normally classified into three sub-bands based on their frequencies, namely low, middle and high frequency-bands. What is the threshold value we can use to classify a DCT Coefficient in low, middle or high frequency band? In other words, how can we classify the DCT coefficient frequency bands radially? Below is an example of the radial classification of the DCT coefficient frequency bands.

The idea is based from the paper : Noise Characterization in Ancient Document Images Based on DCT Coefficient Distribution


texts:
I'm attempting to design a FIR high pass filter than keeps signals above 200Hz and rejects signals below 60Hz with a sampling frequency of 500 samples/sec. This is my first time attempting this and I'm a little confused.

I started with looking at pole/zero placement. I'm not exactly sure how to go about this but this is what I did. I want a zero on the unit circle at the 60Hz point to minimize that frequency. To find this spot on the axis I did the following:

$$\frac{60}{\frac{500}{2}}\pi=0.24\pi$$

So I believe I need a zero at 0.24pi on the unit circle on the z-plane. Does this sound right? How would I then transform it in a difference equation in the sequence domain?


texts: Looking for pratical quantitative comparison metrics for scaled, delayed and warped Signals

texts:
[Nota: I have validated the only answer so far, yet more explicit versions are still welcome]
The following question is detailed in 1D, with time as the ordinal variable. Similar questions could apply in other dimensions.
In several signal processing techniques, such as blind source separation (BSS), filter banks, or deconvolution, one may wish to estimate a signal $x(t)$ and only recovers $s.x(t+d)$, a scaled and delayed estimate. Rotations and shears can be added in higher dimensions, and many other more.  $s$ is a scale factor, $d$ a delay. One might even stumble with warped data ($x_{s,d,w} = s.x(t/w+d)$), as  in super-resolution for instance.
In theory, one can estimate continuously $s$ and $d$ with local correlation or Fourier transforms (How to match 2 signals that have same information, though shifted and scaled). The warping $w$ might be estimated with the scale-transform or wavelet representations. I have read several BSS papers and books, asked people, been in conferences, and was not able to find a standard, or at least a usable metric.

In image (it works on signals as well), the Structural Similarity index somehow compensates offset and variance.

Are there practical error metrics to compare the original $x(t)$ with the transformed $x_{s,d,w}(t)$ in the context of sampled signals and noise conditions? Indeed, the discretization induced by sampling complicates the comparison task (imagine for instance a $1$-sample spike on the sampling grid that would be delayed by a non-integer time), as well as the noise.
Should one resort to asymmetric quantities such as divergences?
Can other signal properties help (bandpass, sparse, positive, etc.)?

Forgetting about the warping, I have tried to minimize a standard $\ell_p$ norm, with $s$, $d$, and $w$ as parameters, and to smooth both signals. I am not satisfied with the complexity and the outcomes, and this is a bit tedious.


texts: Correcting large variable dc bias in signal FFT

texts:
I am somewhat new to signal processing, but i'm hoping I can provide enough information in order to be pointed in the right direction.

I am working on analyzing the frequency spectrum of audio signal impulse response.  A few of the impulses seem to have what looks like a large variable dc bias, and when performing a FFT they come back with very large DC values.

subtracting the mean did not help much in ways of correction.  I think there may be other issues besides what appears to be DC but I am stuck.

The graphs below are plots of the FFT results from two signals - one with the warped dc bias, and one without. Hopefully this demonstrates better what may be happening to cause this.  Any help is greatly appreciated.

both Signals are of the same sample length and FFT size

1st Image - FFT of good Signal, small-no DC value

2nd Image - FFT of problem Signal Large DC value


texts: What is meant by Beamspace domain and antenna domain for channel matrix?

texts:
I am not from Engineering background that's why asking these overly simplistic question. Anyone who can answer. These are basic research questions.
Thanks


texts:
I made a terrible vocoder earlier and got some help on it here:

My vocoder sounds bad

My new implementation is much better. I'm using the band pass filter described in the document suggested on the old question. It seems to work much better, but the carrier wave is bleeding heavily into the output signal to the point where I can faintly hear some good vocoder sound effects behind the sound of the carrier wave. This bleed gets louder when the input volume is higher. Is this a common and easily fixed problem or should I give more info for you guys to help me debug this?


texts: Improving Template Matching Algorithm for OpenCV

texts:
Car detection using OpenCV
I am working on a small personal project where i have to know wether an image shown is a car or not. I am using the basic OpenCV python based template matching.
Positive Image / Template Image
Car Top View  :-

The simple template matching by using one of the positive image on the other is giving the required result.
Negative Image
But when we are using negative images like

the template finder is finding this as a positive match.
Methods tried and failed :-

Increase the threshold for template matching :- Increasing the threshold is causing many of the actual positive images to stop matching

Different types of template matching techniques that are there in OpenCV :- Not giving any better result

Using multiple templates to reinforce positive and negative matches:- Doesn't work well for all the cases

Is there  a better way to remove the template matching false positives that we are getting. Are there any feature matching or edge detection based techniques that can be used instead of template matching to improve my algorithm.


texts: How to evaluate an impulse/response model in MATLAB

texts:
This is my first experience in signal processing, so don't treat me badly :)
However, I have estimated an impulse-response model using the MATLAB System Identification toolbox. The output, validated on unseen data, has a good fit, approximately around 80%. However, when I export the object imp - its type is idtf - containing the estimated model, I don't know how to evaluate it on unseen data. Obviously, in the structure is contained the polynomial structure of the model, but it bothers me to write a function to estimate it... Is there any quicker procedure?


texts: Identifying the correct source of packets

texts:
I have a basic question regarding digital communication , especially pertaining to wireless internet .

A packet is basically a stream of bits . And a bit is detected as zero or 1 depending on whether the amplitude of a signal of certain frequency is crossing a threshold or not .
Now suppose a particular host or source $ A $ is sending this signals at a frequency $\ f $  to a receiver $ R $
Then it is possible that some device , $ B$ near the receiver can be placed to detect the frequency and regenerate a signal in that frequency and direct it towards the receiver .

In internet communication , I guess the packets will definitely be containing information about the source IP address, the IP of $ A $  but $ B $ can alter the bits so as to insert send its own IP  address in the packets .

In such situations , how do real receivers , identify the correct source ?


texts: difference beetween decimate an array and decimate function

texts:
i'am actually comparing data and i got a question that i cannot anwser,
here is what i understand.
nanometrique move
X[::256] : decimate and take 1/256 point in the vector (blue)
decimate(X, 256, ftype='fir') : decimate and apply a LP filter (orange)
questions are : why those peaks on decimate methode ? and as i know of my lessons a LP filter does not permit a 30dB gap beetween those figure, an idea why ?
for the PSD i apply mlab.psd function

any sugestions ?
up after comment :
thanks for your reply
i'll try to explain clearly,
blue signal is raw data with decimation of 256 without any filter
orange signal is apply with the function decimate on the same raw data and that function decimate apply a FIR
for the LP i don't even know is frequencial response du to the function
the original signal is like 200 000 000 data and it is too heavy for my computer ...i reduce the raw decimation : 256 to 29


texts:
I'm new to signal processing and I'm learning on my own. Therefore I might ask questions not worthwhile your time.

At the moment I'm trying to figure out the concept of a lowpass filter.
Using matlab I created a complex signal x:

x = rand(1,1000)+rand(1,1000) .* 1i

And i created a filter Num with the fdatool, the specifications of Num are :

Fs = 1000;  % Sampling Frequency

Fpass = 460;             % Passband Frequency
Fstop = 500;             % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor

From what I understood I had to do now. I tried to do the Convolution of Num and x in the frequency domain by multiplication of the vectors.  Step by Step I did this:

First I padded Num to the same length that x has.
Then I multiplied x and Num, to get y= Num * x;

Now y should have nearly no frequencies above 500, right?
If I check the signal there doesn't seem have happened a lowpass filtering.


texts: Why is Fourier space not adequate for (theoretical or digital) filters?

texts:
As far as I have seen, almost all theoretical filter design occurs in Laplace or Z-space. Also, there is a pervasive connection to real life analog filters in the design. If one is just thinking in a mathematical theoretical thing (or something that could be implemented digitally), why wouldn't one filter signals in Fourier Space?
Why is, say, multiplying the Fourier Transform of a certain function by a unit step up - step down function, and then making the Inverse Transform of the resulted signal a "band pass filter"? Why should one use Butterworth filter or similar things to make a digital filter?


texts: Difference between these two Continuous Wavelet Transforms?

texts:
I am porting Synchrosqueezing to Python, and came across an implementation difference on CWT between mine and PyWavelets' - details below. The idea is to merge this implementation to PyWavelets if latter doesn't already have it, but I cannot explain the difference as I'm unfamiliar with wavelet transforms.

What is each implementation trying to accomplish, and what are some key differences?

Full post; comparing pywt.cwt vs. cwt_fwd, pasting transform excerpts:

# pywt.cwt
for i, scale in enumerate(scales):
    step = x[1] - x[0]
    j = np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step)
    j = j.astype(int)  # floor
    if j[-1] >= int_psi.size:
        j = np.extract(j < int_psi.size, j)
    int_psi_scale = int_psi[j][::-1]

    size_scale = next_fast_len(
        data.shape[-1] + int_psi_scale.size - 1
    )
    if size_scale != size_scale0:
        # Must recompute fft_data when the padding size changes.
        fft_data = fftmodule.fft(data, size_scale, axis=-1)
    size_scale0 = size_scale
    fft_wav = fftmodule.fft(int_psi_scale, size_scale, axis=-1)
    conv = fftmodule.ifft(fft_wav * fft_data, axis=-1)
    conv = conv[..., :data.shape[-1] + int_psi_scale.size - 1]

# cwt_fwd
xh = np.fft.fft(x)  # x == data
xi = np.zeros((1, N))
psihfn = wfiltfn(wavelet_type, opts)

for ai in range(na):
    a = Wx_scales[ai]
    psih = psihfn(a * xi) * np.sqrt(a) / np.sqrt(2 *PI) * (-1)**k
    dpsih = (1j * xi / opts['dt']) * psih

    xcpsi = np.fft.ifftshift(np.fft.ifft(psih * xh))
    Wx[ai] = xcpsi

    dxcpsi = np.fft.ifftshift(np.fft.ifft(dpsih * xh))
    dWx[ai] = dxcpsi

Observations:

cwt_fwd has an imaginary component, pywt.cwt doesn't
pywt.cwt's norm is ~10x greater (see norm_abs) [in 'full post']
The signal (x) preprocessing seems to undo itself in this case as evident in vs. Unprocessed
pywt.cwt handles the noisy signal (sN) case considerably differently in the upper-half

Visualized differences:


texts: Filter design, relationship between energy and stop-band ripple

texts:
I read some example of design LPF which I didn't understand something.
The stop-band in that example is $\frac { 22 }{ 25 } $  from the over-all frequency, and I want to filter some white noise.
After we found out in that example that the energy of the white noise in the stop band area is $0.22\%$ from the overall noise before filtering, they used the equation:
$$\frac { 22 }{ 25 } { \delta  }_{ s }^{ 2 }=0.0022$$

Why ${ \delta  }_{ s }^{ 2 }$?
How does is related to power/energy relationship?

The full example in the book of Boaz Porat page 247, can I upload a picture from a book?


texts:
I'm trying to make some simply autofocus algorithm for my device and pretty confused with 2D FFT results.
As I know for now - number of high-frequency's components in 2D FFT will say me that image is more focused, than other one where 2D FFT got more low-frequency's components. And i do not know how to calculate this parameter.

I'm thinking about applying high-pass filter, but how i should choose cut-off frequency? In 1D FFT it's pretty simple to determine just looking at the picture of 1D FFT..


texts:
Are there general techniques to derive DTFTs? Given a bandlimited function $x(t)$, how do I find

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

Generally, it is easier to derive the continuous transform (never mind the constants):

$$X(f)=\int_{-\infty}^{\infty}x(t)e^{-i \omega t}\, dt$$

because we have a wealth of integration theory to fall back on. Ok, one can also calculate $X(\omega)$ from $X(f)$ by adding shifted copies, but we're again back to summing an infinite sequence. Any ideas on how to approach DTFT derivations in a smarter way? Or is sweating out the tiring summation the only way? I don't have any particular function in mind, but if one has to be chosen as an example, I'd pick $x(t) = \text{sech}(a t)$. The continuous FT of this is known and can be found in the table on Wikipedia (#208).


texts:
I am getting started with Python's FFT. I tested it on a signal that is a sum of three signals, two of which have an eigenfrequency of the grid, the third one does not (but due to large no. of data points, I still get a peak at the correct position).

However, the amplitudes are a bit off. The amplitude of the third signal is a bit higher than it should be (>=5.1 instead of 5). I was expecting that the amplitude would be lower for this signal, as it does not have one of the eigenfrequencies. But how can it be higher?

Also, for the two "well-defined" components of my signal, the amplitude is only correct to order 10^(-2). Again, one of them is a bit larger than the actual value. Are these numerical errors in the fft routine or is something wrong with my code?

Thank you!

import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fft, fftfreq, ifft

L = 5
N = 1000
pi2 = 2*np.pi

x = np.linspace(0, L, N)

y1 = 10*np.sin(6*(2*np.pi/5)*x)  ## has one of the eigenfrequencies
y2 = 20*np.sin(10*(2*np.pi/5)*x) ## has one of the eigenfrequencies
y3 = 5 *np.cos(626.3*x)  ## lies between the highest and second highest eigenfrequencies

y = y1 + y2 + y3

A = fft(y,N)

freqs = fftfreq(N, d=L/N) * pi2
mask = freqs > 0

plt.figure(1)
plt.plot(freqs[mask], 2*np.abs(A[mask]/N))
plt.show()


texts: Does modified Morlet wavelet function satisfy the admissibility condition?

texts:
Does modified Morlet wavelet function satisfy the admissibility condition ?
Do you have any reference for the answer please.
Thanks.


texts: Colored noises: Stationary or non-stationary?

texts:
I know for sure that white noise is considered as a stationary sound, but is that true for the rest of them? http://en.wikipedia.org/wiki/Colors_of_noise
Also, can we meet some of these sounds in nature? Is there any relation between these noises and noises that you can hear in a town(vehices,traffic,people chatting) or are these noises just an artificial invention?

Thanks in advance for your ideas and thoughts!:)


texts: What's the suitable specifications for the FIR filter?

texts:
i have EEG signal,the sampling frequency is 256 HZ .The wavelet decomposition that i will use does not allow the extraction of specific frequency bands without additional filtering. Thus, in order to correlate the wavelet decomposition with the five sub-bands of the EEG signal,the frequency content segments of the EEG signal was restricted to the 0-64 Hz band by convolving the signal with a low-pass finite impulse response (FIR) filter.So when i try to design the FIR filter using MATLAB i want to know how to choose the suitable specifications.

thanks


texts: Frequency of Signal of a point tracked in a video?

texts: I have used MATLAB to track a feature point (corner point) using Lucas-Kanade tracker and plotted its movement through each frame now I need to find its frequency content, how should I find it? (I have only taken $y$-axis movement of the corner point)


texts: A basic question regarding frequency analysis of an EEG signal

texts:
Assume that an EEG signal is sampled at $f_s = 300$ Hz then a 10000-point segment of it is selected, called $x[n]$. The corresponding 10000-point DFT is then computed and called $X[k]$.
Assume further that famous EEG frequency bands are known as the following intervals (in Hz):
$ \delta=[1,4],$ $ \theta = [4,7],$ $\alpha = [7,12],$ $ \beta = [12,30]$.
Please help me decide which of the following 4 statements are true, and which are false?

Any 5000-point of X[k] is enough to obtain the entire frequency information.

My attempt: According to the famed Nyquist criterion, the maximum frequency we can have is $f_m = f_s / 2 = 150$Hz. The corresponding index $k$ can be found this way:
$$ \frac{k f_s}{N} = f_m = \frac{f_s}{2} \rightarrow k= \frac{N}{2} = \frac{10000}{2} = 5000.$$
This implies that the maximum frequency happens at the 5000'th frequency index. But it doesn't mean any 5000-point segment gives the frequency information so this statement is incorrect.

$X[k]$ for 1000<k<1500 is enough to get the $\alpha$ band information.

My attempt: At $k=1000$ we get
$$ f = \frac{  k f_s }{  N } = \frac{  1000 \times 300 }{  10000 } = 30 \textit{Hz}$$
and at $k= 1500$ we have
$$ f=\frac{ 1500\times 300  }{10000   } = 45 \textit{Hz}$$
which are out of the $\alpha$ band range, so this statement is incorrect as well.

$X[k]$ for 8000<k<10000 is enough to get the $\alpha,\beta,\delta, \theta$ bands information.

My attempt: For the very same seasons just discussed, $k=8000$ corresponds to $f= 240$ Hz and $k= 10000$ corresponds to $f= 300$ Hz  which are out of all bands' ranges. As a matter of fact, 240 and 300 Hz are above the maximum frequency $f_m$ allowed by Nyquist's criterion, so I don't really know how to interpret this. Maybe I'm mistaken?

A lowpass filter with frequency $ \frac{\pi}{2}$ removes power line frequency at $60$ Hz.

My attempt: I think this is trivially correct because $\frac{\pi}{2}$ Hz is about 1.57 Hz and the lowpass filter zeros out anything above this frequency, including 50 Hz.


texts: Cochleagram vs STFT for CNN-based speech segregation

texts:
I’m working on a project (python-based) that would use ideal ratio masks (IRMs) as a basis for cleaning noisy speech in various environments. Specifically, this will be accomplished through the use of a convolution neural network, where input features would be either a cochleagram or STFT data matched with its respective IRM label.

(This is the article I’m basing my project on --> http://web.cse.ohio-state.edu/~wang.77/papers/CWYWH.jasa16.pdf)

As I’m (relatively) fairly new to learning about digital signal processing, I was hoping someone would be able to give me an explanation hilighting some of the differences between cochleagrams and STFTs, and what their individual use-cases would be.

Many thanks in advance!


texts:
I computed spectrum analysis of sine wave (without windowing and zero-padding) and I got plots as follows:

As far as I know they are correct. When I learnt about DSP, I got the huge problem with understanding what is the phase spectrum. Recently I realised that "phases" are translations along time axis of each sinusoidal component. So I don't see any point to get them out of range $\langle -\pi, \pi \rangle$. I computed my spectrums with python and to get the phase spectrum I used unwrap function from numpy. As I understand this documentation correctly, it is impossible to get the result larger than $\pi$ or less than $-\pi$ unless I set the discount parameter. It does make sense for me.

Could you explain me:

whether I understand the phase spectrum correctly?
why the values of phase spectrum are out of range $\langle -\pi, \pi \rangle$?
how does the function unwrap work?


texts: How to encode a byte array to sound file and then convert back?

texts:
This question is moved from Stack Overflow.

Sorry about this noob question, because I never work with matlab and signal processing before. Here is what I want to do: I have a fixed length of byte array X, now I want to encode it to a sound file, I also want this process to be reversible, which means the sound can be converted back to X with no error. I searched online, and found the following code:

M = 16;
x = randint(5000,1,M);
y=modulate(modem.qammod(M),x);

My question is that, is QAM the best way to do this? and how to use it? I will need to send the sound signal over the air, and I just need to transfer about 128 bit data each time. A little bit code example will be really appreciated, Thank you!

update#1: I tried to output y by sound(y), but matlab does not allow me to do so, it says I can only output floating numbers. How can I solve this? Thank you!


texts: Get data after "Correlate Access Code Tag"

texts:
in my Gnuradio flowgraph i have a stream of bytes that are 0x00 or 0x01.
Using the "Correlate Access Code Tag" block, I am able to insert a tag just after that a specific Access Code has been found (Say 11111101010101x the tag is on x if Access Code is 11111101010101)
Now i want to insert into a vector only a specific number k of items starting from x.
Say k=8, if the stream is  000100010001111110101010100110001 the vector should contain
00110001.  If i use the "stream to vector" block after the "correlate access code tag" block with vector length=8 I end with a wrong division since the stream to vector seems to operate on the whole stream from the beginning and does not take into account the correlate access code tag block.  In my example the vectors i have resulted divided in this way  00010001 00011111 10101010 10011000 1
So what is the right way to pick n items after correlate block ??
Tnx, Fabio.


texts: Last year Apple releases their new phone Apple X that has an exciting feature called Animoji. I am so fascinated by this technology that I really want to make my own version of it. But the problem is I am really new in image processing stuff and I've no idea about the technology behind Animoji. So what is the technology behind the Applle Animoji and what sort of knowledge should I have if I want to build my own version of Animoji from scratch?


texts: What is the right way to calculate the inverse Z-transform of $zX(z^{-1})$

texts:
say the signal $x(n)$ has the z transform $X(z)$ and there is signal $x_1(n)$ that

$X_1(z)=zX(z^{-1})$

I tried 2 different approach to get the relationship between $x(n)$ and $x_1(n)$ and the result is different without any error being noticed

Approach 1:

$x_i(n)=x(-n) \rightarrow X_i(z) = X(z^{-1})$

$x_1(n)=x_i(n+1) = x(-n-1) \rightarrow X_1(z) = zX_i(z)=zX(z^{-1})$

Approach 2:

$x_i(n)=x(n-1) \rightarrow X_i(z) = z^{-1}X(z)$

$x_1(n)=x_i(-n)=x(1-n) \rightarrow X_1(z) = X_i(z^{-1}) = zX(z^{-1})$

Hence in the end, I don't know $x_1(n) = x(-n-1)$ or $x_1(n) =x(1-n)$ and I can't find any mistake. I know that different function can have same z-transform but different ROC but in this case, no any info of ROC is given


texts: How do I efficiently evaluate a convolution integral between different-sized domains?

texts:
This is a question that I've previously asked over on math.stackexchange, and I have yet to receive a useful answer. It was suggested that I post this here.

The problem itself originally comes from simulation of wave propagation, and involves large two-dimensional arrays of complex data. I'll explain it here as a simplified one-dimensional problem, however.

Suppose I have the following integral relationship involving a convolution:
$$
f(x) = \int_{-\infty}^{+\infty}dx'\, g(x')\, K(x - x')
$$
Assume that the functional form of the kernel $K$ is known, and suppose that the function $g$ has compact support on $[-a/2, +a/2]$ for some width $a>0$. Assume further that I have data about the function $g$ in the form of $N$ equally-spaced samples $g_n$:
$$
g_n = g(x_n)\, ,
$$
where
$$
x_n = -\frac{a}{2} + n\frac{a}{N}\, , \qquad n = 0,..., N-1
$$
My goal is to generate some $N$ evenly-spaced samples of the function $f$.

If we suppose that we want to evaluate $f$ on the same domain on which $g$ lives, and for which we have samples of $g$, then it's easy. I discretize everything in the usual way, and the whole thing becomes a discrete convolution:
$$
f_n = f(x_n) = \sum_{m = 0}^{N-1} g_m K_{n - m}\, ,
$$
where
$$
K_{n - m} = K(x_n - x_m) = K\left(\frac{a}{N}(n-m)\right)\, .
$$
Everything is simple since I can use the discrete convolution theorem: I FFT the $\{g_n\}$ and $\{K_n\}$ sequences, multiply them, do an inverse FFT, and I'm done. I have my $\{f_n\}$ data. The whole thing will be an $O(N\log N)$ process.

Suppose instead that I have to evaluate $f$ on a different-sized domain. For example, perhaps the $g$ data is narrow, but the kernel $K$ is very wide. Let's say I have to generate $N$ samples of $f$ on the domain $[-L/2, +L/2]$, for some domain width $L > a$. $L$ could even be much larger than a.

If I try to discretize everything in a similar way, I'd have to first define a sequence of points $\{X_n\}$ in the wider domain:
$$
X_n = -\frac{L}{2} + n\frac{L}{N}\, , \qquad n = 0,..., N-1
$$
Then the data I seek is $\{f_n\}$ where $f_n = f(X_n)$. The original integral is still over the smaller domain, so we have:
\begin{align}
f_n
&= \int_{-a/2}^{+a/2}dx'\, g(x')\, K(X_n - x') \\
&\approx \sum_{m = 0}^{N-1} g_m\, K(X_n - x_m) \qquad \text{(discretize)}\\
&= \sum_{m = 0}^{N-1} g_m\, K\left(\frac{L-a}{2} + \frac{a}{N}(\alpha\, n - m)\right)\, ,
\end{align}
where $\alpha = L/a > 1$ is the ratio of the two domain sizes.

This is no longer a discrete convolution, since the second term is no longer a function of $n-m$, but instead of $\alpha\, n - m$. So I can no longer use the convolution theorem. I can imagine evaluating the sum above by brute force, but that would be an $O(N^2)$ process, which is prohibitive for my problem.

Is there any way of efficiently calculating the discrete version of this convolution integral between two different domain sizes? Thanks.


texts: GNU Radio head block doesn't finish simulation

texts: My task is to generate finite number of samples and save it to the file. For that purpose I used head block but after setting the run option to: Run to Completion, the flowgraph doesn't stop, just hanging. Any idea how to force my flowgraph to stop after certain amount of samples is written to the file?


texts: Influence of image size to edge detection in OpenCV

texts:
I need to scan an image in OpenCV and obtain the lines. For this I use the canny edge detector and then the propabilistic Hough transform, with which I actually get the lines.
The problem now is that I get many short lines instead of a few continuous ones. The image I am analyzing has the size 1024x768 and I am wondering if this could be an error source. Theoretically making the image smaller should lead to fewer lines. Might that help me in my task? Or are there any other factors which I should take into account when extracting lines?


texts: Tied rank adjustment for Spearman Correlation

texts:
I have two 16 sample vectors, a and b as listed below:

a = 345.86  349.83  257.84  438.43  195.03  0   341.39  0   0   0   0   0   0   0   0   0

b = 313.30  349.83  0   394.07  195.03  257.84  0   0   0   0   0   341.39  0   0   0   0

My aim is to calculate the Spearman correlation between the two vectors. I know that it can be done easily in MATLAB by using the function:

corr(a,b,'Type','Spearman');

However, my aim is to calculate this step by step, i.e.

Calculating the ranks of vectors a and  b
Applying the formula $C=1−\frac{6\Sigma(R_x−R_y)^2}{n(n^2−1)}$
Adjusting the tied ranks (most important part)

Firstly, using the function corr(a,b,'Type','Spearman') I get the correlation coefficient as 0.5845.

Now as to calculate the same step by step, I am using the tiedrank function to calculate the ranks for the vectors a and b. From The Mathworks documentation:

  [R,TIEADJ] = tiedrank(X) computes the ranks of the values in the vector X. If any X values are tied, tiedrank computes their average rank. The return value TIEADJ is an adjustment for ties required by the nonparametric tests signrank and ranksum, and for the computation of Spearman's rank correlation.

Using this for the above example, I get the following:

Rank_a = 14 15  12  16  11  5.5 13  5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5

Rank_b = 13 15  5.5 16  11  12  5.5 5.5 5.5 5.5 5.5 14  5.5 5.5 5.5 5.5

And the Tied Adjustments that I get:

TIEDADJ_a = 495

TIEDADJ_b = 495

So, later I calculate the differences of squares of the ranks between a and b i.e. Sum of (Rank_a(i) - Rank_b(i))^2 which is equal to 214.

I do not know how to take into account the TIEDADJ_a and TIEDADJ_b values in addition to the above summation of difference of ranks calculated, so that I can achieve the correlation coefficient calculated by the MATLAB function directly.


texts: Confused with using Fourier Transform Properties

texts:
I'm seeking the Fourier transform of $t \cdot (u(t+1)-u(t-1))$

Given the fact that (The first minus here should be a plus)

Using properties.
My initial approach was:

but I wasn't sure how I should go around nesting properties, so I decided that I can perhaps:

However, checking Wolfram, we have:

I don't understand what's going on, my questions are:
1- What went wrong in my second attempt?
2- How would I carry on if I chose to go with my first attempt?
Sorry if these questions seem naive or anything, I've just been so recently introduced to the topic.


texts: YIN Pitch detection Algortithm ( how do I improve my results )

texts:
I am using YIN algorithm in a school project of mine which uses pitch detection on guitar sound. I when I play a note I get random frequencies at the beginning until they stabilize. I am thinking those are probably from action of pick on the strings.

I am going through the original paper:

  Cheveigne A, Kawahara H. - YIN, a fundamental frequency estimator
  for speech and
  music

trying to reverse engineer the library and improve my results.
I am studying computer science and my knowledge of signal processing is limited.
A summary:

Step 1)
Auto-correlation :- We try to find the correlation of the signal with itself shifted by a lag within a window.

Function can possibly have infinite values. We chose the highest peak with non-zero lag. Within a lag range.(Why the highest peak? Does it mean the loudest frequency?). The paper says if upper limit for $\tau$ is high. The algorithm may chose higher order peak (what are higher order peak?)

The following steps are to improve the accuracy of the results

Step 2)
Difference Function :- Model the signal in form of a difference function.
$$ d(\tau) = \sum_{j=1}^{j=W} (x_j - x_{j+\tau})^2 $$. Which gives :
$$ d(\tau) = \sum_{j=1}^{j=W} (r_t(0) + r_{t+ \tau}(0) -2r_t(\tau)$$

So we're basically using amplitude as bias.

step 3)
Cumulative mean normalization  Replace the difference function by cumulative normalized difference to avoid selecting value with zero lag:

$$d_t(\tau')=\frac{d_t(\tau)}{(1/\tau)\sum_{j=1}^{\tau}d_t(j)}$$

step 4)
Absolute Threshold ( Could Anyone explain this section?)

step 5)
Parabolic Interpolation : Fit the $d(\tau)$ estimates to a parabolic curve.

step 6)
Choose The Best Local Estimate : Self explanatory

I am trying to compare guitar sound with a monophonic midi.

I think the parameters I should be thinking about tinkering with are window size and threshold to improve my results or I could discard first few frames.
Could anyone point me in the right direction?

The parameters I am using :

SAMPLE RATE : 44100

WINDOW SIZE : 1024

HOP SIZE : 512

THRESHOLD : 0.1


texts: How to simulate (in matlab or c) power spectrum of random brownian noise signal?

texts:
I have a semi random noise signal, it does have some periodicity however.

I can estimate the power spectrum (of the blue signal, I thought that it would be more clear, because white noise is flat, but I digress)

How could I generate a noise signal from a power spectrum?

Brownian noise estimation is not the same as estimating from a power spectrum because of the 'sine wave' noise in the spectra, I would like to capture that.


texts: How to find the solution of nonlinear least square (NLS) equation for freqeuncy estimation

texts:
I am struggling to understand the theory of parametric methods for power spectral estimations of line spectra.

And I want to find the solution of nonlinear least square (NLS) method for estimating the frequencies of line spectra. For example,

$Let,$
$$f(\omega, \alpha, \phi) = \sum_{t=1}^{N} |y(t) - \sum_{k=1}^{n} \alpha_{k} e^{i(\omega_{k}t+\phi_{k})}|^2,$$

$then,$

$$f = (Y-B\beta)^*(Y-B\beta) \\= [\beta-(B^*B)^{-1}B^*Y]^*[B^*B][\beta (B^*B)^{-1}B^*Y]+Y^*Y-Y^*B(B^*B)^{-1}B^*Y $$

$where$
$$Y = [y(1) ~y(2) ~... ~y(N)]^T$$
$$ \beta = [\beta_{1}~...~\beta_{n}]^T$$
$$B =\begin{bmatrix} e^{i\omega_{1}} & ... & e^{i\omega_{n}} \\ \vdots &  & \vdots \\ e^{iN\omega_{1}} & ... & e^{iN\omega_{n}} \\ \end{bmatrix}.$$

The best choice of the $\omega = [\omega_{1}, \cdots, \omega_{n}]^T$ to minimize $f$ is,

$$\hat \omega = \underset{\omega}{\operatorname{argmax}} [Y^*B(B^*B)^{-1}B^*Y]$$

I have two questions for the above equations.

1)How can we derive $$f = [\beta-(B^*B)^{-1}B^*Y]^*[B^*B][\beta (B^*B)^{-1}B^*Y]+Y^*Y-Y^*B(B^*B)^{-1}B^*Y ?$$
2)How can we solve the $\operatorname{argmax}$ problem for $n$ different parameters of $\omega$?