texts: Effects of averaged periodogram (Welch's Method) on resolution?

texts:
I'm trying to make an app(running on an android phone) that can detect the frequency of a guitar string being played. I know it'll be hard to detect that but that's a problem I'll deal with later.

Right now I'm trying to reduce the effects of noise on my FFT of an audio signal. I'm sampling at $44.1\textrm{ kHz}$ and have a sample size of $4096$. This gives me a bin size of $10.7\textrm{ Hz}$.

I've read good things about averaged periodograms killing noise. I want to try to use an averaged periodogram, Welch's Method, to reduce the effects of noise.

What I'm wondering is, if I window it into $1024$ sample windows, then will the resolution of my final product become $\approx 40\textrm{ Hz}$?
If I have overlapping windows to keep my resolution high, will I still have effective noise reduction?

Sorry if this is a bad question. This is my first time doing DSP outside MATLAB so I'm struggling.


texts: Calculating energy spectral density function of half-cosine function

texts:
I have the following function:

$$ g(t) = A\cdot \textrm{rect}\left(\frac{t}{T}\right)\cos\left(\frac{\pi  t}{T}\right)$$

How do I find the energy spectral density function? I think? it is defined as :

$$ \int\limits_{-\infty}^\infty \left|g(t)\right|^2\, dt $$

After squaring it, I'm a bit lost since I do not know how to continue...


texts: What is different between gradient and imgradient functions in MATLAB

texts:
I would like to compute the gradient of an image I. I have the following two options :

[Gx, Gy] = gradient(I);
g = sqrt(Gx.^2+Gy.^2);

or

[g,~] = imgradient(I, 'sobel');

My questions are

What is the gradient method used in the first option ?
What are the benefits of using Sobel method to estimate the gradient ?


texts:
I'm trying to reproduce this !!!WARNING LOUDNESS!!! sound in C. In the website, it says to use Carrier: 500 Hz, modulator frequency: 1 Hz, modulator index: 500, however I'm getting lots of clicks and fuzzing (along with the right sound) with this code:

#include <stdio.h>
#include <math.h>
#include <string.h>

#define DURATION    20
#define RATE        44100
#define NUMCHANS 1

#define NUMSAMPLES  (DURATION*RATE)
#define NUMSAMPLESC (NUMSAMPLES*NUMCHANS)

#define BYTES_IN_HEADER (11*sizeof(int))
#define BYTES_IN_BUFFER (NUMSAMPLESC*sizeof(short) + BYTES_IN_HEADER)

#define DELTAT (1.0/RATE)

#define M_PI 3.14159265359

static const int wavHeader[] = {
0x46464952,
NUMSAMPLESC*2+36,
0x45564157,
0x20746D66,
16,
1|(NUMCHANS<<16),
RATE,
RATE*NUMCHANS*sizeof(short),
(NUMCHANS*sizeof(short))|((8*sizeof(short))<<16),
0x61746164,
NUMSAMPLES*sizeof(short)
};

static short buffer[BYTES_IN_BUFFER/sizeof(short)];

int main(void) {

 unsigned int i;
 double t = 0.0;
 FILE *fp;

 memcpy(buffer, wavHeader, BYTES_IN_HEADER);

for(i = BYTES_IN_HEADER/sizeof(short); i < BYTES_IN_BUFFER/sizeof(short); i++, t += DELTAT) {
    buffer[i] = 100 * sin(2*M_PI*500*t + 500*sin(t*2*M_PI*1));
}

 fp = fopen("test.wav", "w");
 fwrite(buffer, 1, BYTES_IN_BUFFER, fp);

fclose(fp);

return 0;
}

I have tried using single precision floats, reversing endianness, but with the same result. The only thing that seems to remove the clicks and fuzzing to a certain extent (but they're still there) is if i copy the lower 8 bits of each sample into the upper 8 bits.

My sound card is integrated from a P43DE3 motherboard and can handle audio CDs well, and my OS is Windows 8.

How do I reproduce that sound without noise?

EDIT

It seems to produce the correct output when I run it on a Mac. Perhaps and endianness error? Will investigate further.


texts: How can you 'perfectly' reconstruct phase for transients in a modified audio absolute spectrogram?

texts:
I am currently working on the absolute spectrogram of an input audio signal. I am using non-negative matrix factorization to separate out the rhythmic components (bass drum, snare drum and hi-hat) of the file and shift around the activations of these components to get a modified spectrogram.

Now, I want to reconstruct the phase for this modified spectrogram so that I can take an inverse FFT and recreate the audio signal in time domain. I am currently using the Griffin Lim and the Phase Gradient Heap Integration methods from this phase retrieval toolbox. Unfortunately, this doesn't produce good results. Maybe it's with the arguments I'm using. I'm currently using a window size of 2048 with a hop of 512.

1, 2 - These questions were previously posted where the idea of reconstruction of audio from absolute spectrogram was tackled. This question pertains to how well one can reconstruct the audio signal back in the time domain and what would be the best way to do so when the harmonic content would be the same but the transients would be time-shifted. Also, how would you objectively or subjectively evaluate such a system formally?

Thanks!


texts: Linear phase Notch filter / Band reject filter implementation in C++

texts:
I'd like to implement a digital notch filter, which has linear phase, in C++ or C#.
Example : a notch filter that removes 400hz frequency on a .wav file (16 bit, 44.1 Khz, stereo).
I'm a bit lost with many examples, I don't achieve to implement them.
Maybe some of you has a working example?

Thank you in advance.


texts: finding inverse Z transform with usage of tables (LTI and causal sequences)

texts:
problem is as follows

try to find the inverse Z

$Z^{-1} (\frac{3}{z+2}) = ???$

with the usage of z-transform tables

Ok, so in order to find something from the table, I thought that we expand with z and try to see if we  get closer to any kind of expression in the tables

$Z^{-1} (\frac{3z}{z*(z+2)}) = Z^{-1} (3*\frac{1}{z}\frac{z}{z+2}) $

then we get constant out

$=3*Z^{-1}(z^{-1}*\frac{z}{z+2})$

then we see that right side fraction can be found in the table, but the combined sum expression can also be found with further manipulation (as can be seen later)

So, I think we declare temp variable here to clean it up (this portion is found in the table green portion)

$F(z) = \frac{z}{z-(-2)} \iff f[k]=(-2)^k $

then we can make it so

$3*Z^{-1}(z^{-1}*F(z))$

then we find the direct expression in the table (we have delay of one unit because of $z^{-1}$ term, it is the orange portion)

$3*Z^{-1}(z^{-1}*F(z)) = 3* u[n-1]*f[n-1]$

so it looks like answer should be

$3* u[n-1]*f[n-1] = 3*u[n-1]*(-2)^{n-1}$

my teacher hasnt provided the solutions guide yet, despite that he promised to do so, and I reminded him about it, so I was just checking if this is the right solution and thinking...

I did the same problem, but I used polynomial long-division style to get the following series representation

$3*z^{-1}-6*z^{-2}+12*z^{-3}-24*z^{-4}...$

so based on that I guess that the $x[n]=3*(-1)^n*2^n $

where $ n>=0$

so, I was just a little bit confused does there actually exist the delay of one unit in the original input sequence (result of inverse Z), sometimes I get confused about the causal sequence when you multiply with the unit step sequence, but then, it seems that sometimes also it can just be stated that the original sequence must have been causal sequence, so you just define n>=0 or something like that???


texts:
Could somebody tell me what is the DFT scaling of a $N$-length signal which is padded with $M$ zeros ? Is it $1/(M+N)$ or $1/N$, i.e.
$$
Y_k=\frac{1}{N}\sum_{n=0}^{N-1}y(n)e^{\frac{2\pi nk}{(M+N)}} \quad \text{or}\quad
Y_k=\frac{1}{(M+N)}\sum_{n=0}^{N-1}y(n)e^{\frac{2\pi nk}{(M+N)}}?
$$


texts: The Zero-Crossing rate threshold for a voiced/unvoiced decision

texts: I've implemented a function that calculates the zero-crossing rate for a given signal. I've used this same function to calculate the pitch. To differentiate voiced signals from unvoiced with reference to ZCR: A high ZCR means that the signal is unvoiced and a low ZCR means that it is voiced. My question is whether there is a threshold above which we can consider that a signal is unvoiced.


texts: How to reconstruct the signal which is subjected to FIR filtering - Matlab

texts:
I have a signal $x$ which is corrupted by Gaussian noise.
The signal $x$ is convolved with a filter with coefficients $h = [1, -2, 1]$, to generate signal $y$. Now the signal is taken to Discrete Fourier Transform (DFT) domain by FFT.

In MATLAB,

y = conv(x,h);
z = fft(y);

How can I compensate for the modulation caused by $h$ in the DFT domain.
I have attenuated the frequency components of $x$ by the filter $h$.

I know that to invert an FIR filtered signal, IIR filters are required.
Can I do some compensation (multiplying coefficients of $z$) on Fourier domain.


texts:
In Acero's Spoken language processing (p. 309) it says that the cepstrum of a finite length impulse train
$$x[n] = \sum^{M-1}_{k=0} \delta[n-kN]$$
is given by
$$\hat{x}[n] = \sum^{ \infty}_{r=1} \frac{\delta[n-rN]}{r} - \sum^{\infty}_{l=1} \frac{\delta[n-lMN]}{l} $$

from where is clear that the cepstrum should be non-zero only for values of $n$ that are multiples of $N$.

I'm trying to visualise this result with a quick matlab script ($N =5$ and $M =3$):

delta_train = [1 0 0 0 0 1 0 0 0 0 1];
stem(cceps(delta_train));

but it gives by result

that is inconsistent with the above expression because it's non-zero for every $n$.

I was thinking that this may be somehow related to the way matlab implements cceps() with the FFT and the fact that the computed cesptrum is actually an aliased version of $\hat{x}[n]$, nevertheless, it should still be zero-valued for those $n$.


texts: What's a good way to draw signals and systems diagrams as an instructor?

texts: I'm a TA for an undergraduate signals and systems class. As part of my job, I need to prepare some homework problems. My co-instructors use Microsoft PowerPoint to make the diagrams, but I find it very difficult to use, since we need to copy/paste objects and carefully place them in precise spots. I'd rather write code. Is there some LaTeX package specifically for this purpose? Or, is there an alternate way to draw these diagrams? How do instructors typically do it?


texts: Received Signal Strength Calculation for Probabilistic Line of Sight Channel

texts:
Scenario:
Given a UAV which emits a beacon and the beacon is received by a node on the ground. The node calculates the received signal strength. Considering the channel is Probabilistic LoS, the average path loss is given as:
$$L_{a, k}^{av}(t) = \left[P^{LoS}\eta^{LoS} + P^{NLoS}\eta^{NLoS}\right]\left(\lambda_o d_{a, k}\right)^\alpha$$
Where

$\lambda_o = \frac{4\pi f_r}{c}$
the probability of line of sight is given as:
$$P^{LoS}(t) = \frac{1}{1 + A\exp\left(-B\left[\left(\frac{h_a}{d_{a, k}(t)}\right) - A\right]\right)}$$
$h_a$ is the platform height
$A$ and $B$ are environment dependent parameters
$d_{a, k}$is the distance between the UAV and node
$c$ is the speed of light
and $f_r$ is carrier frequency.
$P^{LoS}, \eta^{LoS}, P^{NLoS}, \eta^{NLoS}, \alpha, f_r$
are the probability of line of sight, mean additional loss for line of sight, probability of non-line of sight, mean additional loss for probability of line of sight, path loss exponent and carrier frequency, respectively.
Assume the transmit power is $P_t$ and received power is $P_r$. And the antenna gains at transmit and receive side are, $G_t$ and $G_r$ respectively.

Goal:
My final aim is to get the distance from platform to node, the node receives the beacon from platform and based on the received power has to ascertain its approximate distance to platform.
Question:
What would be relation for received power at the node, considering the given average path loss above from the platform.


texts: Different ways of decomposing signals into sinusoids

texts:
For purposes of generating a spectrogram, the signal can be divided into overlapping chunks, each chunk multiplied by a window function and transformed to frequency domain using Discrete Fourier Transform (DFT).

But isn't DFT just one of infinitely many ways to decompose a finite-length discrete-time signal into a sum of sinusoids? What other methods are there to do such a decomposition? How would this choice of method affect the appearance of the spectrogram?


texts: Protect an IIR filter from being reverse-engineered

texts:
I created a somewhat unique IIR filter and I want to protect the filter from being reverse-engineered

I think you all know it is quite easy to get all the different weights of an IIR by using impulse signals.

The filter is programatically encapsulated in a bigger program, so I can add some tricks to protect the filter coefficients.

But what would be the best method to protect my IIR filter ?


texts: How to generate realistic echo to sounds?

texts:
I guess the most straightforward way of adding an artificial echo to a sound signal is to add its time-delayed version to the original signal. However, I am not sure how realistic that would sound like.

I was reading papers on room geometry estimation from impulse responses (see for example this paper). So, perhaps one should have several echoes instead of just one.

I guess something also happens to a sound signal when it is reflected (kind of filtering).

What would you suggest here?


texts: Deriviation of the "Twiddle Sum" property

texts:
I can't seem to understand how to derive the "twiddle sum" property:

$$\sum_{n=0}^{N-1}W_{N}^{kn}=N \ \delta[k\bmod N] $$
where $$ W_{N} \triangleq e^{\frac{j 2 \pi }{N}} $$
and $$ \delta[n] \triangleq  \begin{cases}
     1 & \text{if } n=0 \\
     0 & \text{otherwise}
     \end{cases}$$

I've tried doing the following:

$$\begin{align}
\sum_{n=0}^{N-1}W_{N}^{kn} & = \sum_{n=0}^{N-1}e^{\frac{j 2 \pi k n}{N}} \\
 & = \sum_{n=0}^{N-1} \left( e^{\frac{j 2 \pi k}{N}} \right)^n \\
 & = \frac{1 - \left( e^{\frac{j 2 \pi k}{N}} \right)^N}{1 - e^{\frac{j 2 \pi k}{N}}} \\
 & \\
 & = \frac{1 - e^{j 2 \pi k}}{1-e^{\frac{j 2 \pi k}{N}}} \\
 & \\
 & = \begin{cases}
    ?? & \text{if } k=mN \quad\quad m\in \mathbb{Z} \\
     0 & \text{otherwise} \quad k \in \mathbb{Z}
     \end{cases} \\
\end{align} $$

It seems I'm missing an $N$ before the delta, have I done a mistake ?


texts: How does bandwidth affect a modulated digital signal?

texts:
I'm currently reading up on how the lower layers of the ISO/OSI model work and I don't quite understand how broadband communication works yet.

Suppose I want to send some data from point A to point B on a broadband medium. I have assigned a 4kHz wide channel starting at 100kHz. Now I pass some digital data in Manchester encoding over to modem A which transmits the data to modem B.
From my understanding, a Manchester Code signal can just be one of three values at a given point in time - thus, it would require only three different frequencies to be modulated onto.

How would changing the width of the frequency band affect the data rate I can expect? Wouldn't it just stay the same?


texts: Quantization and Sampling - putting it all together

texts:
So after I learned this two topic: quantization and sampling, I'm learning the way to look at both of them and try to optimize the split of a given amount of bit B to N and k, where N is the amount of samples and k is the size of the finite set of numbers can represent the values (quantization).
There is one calculation that I'm not totally understand in the next few slides:

Until now it is very clear, including the yellow formula, let's simplify it more a little bit:
$$=\sum _{l=1}^kP_{\psi }\left(\psi _l\right)\cdot \frac{1}{3}\left(2\cdot \frac{1}{2^3}\frac{\left(\psi _H-\psi _L\right)^3}{k^3}\right)=\sum _{l=1}^kP_{\psi }\left(\psi _l\right)\cdot \frac{1}{3}\left(\frac{1}{2^2}\frac{\left(\psi _H-\psi _L\right)^3}{k^3}\right)=\sum _{l=1}^kP_{\psi }\left(\psi _l\right)\cdot \:\frac{1}{12}\frac{\left(\psi _H-\psi _L\right)^3}{k^3}=\frac{1}{12}\sum \:_{l=1}^kP_{\psi \:}\left(\psi \:_l\right)\cdot \:\frac{\left(\psi \:_H-\psi \:_L\right)^3}{k^3}$$
In the next slide, I don't understand where disappeared:$\frac{\left(\psi _H-\psi _L\right)^{ }}{k^{ }}$. Probably it holds that: $\frac{\left(\psi _H-\psi _L\right)}{k}=\left|\left(delta\right)_l\right|$, but I don't see clearly why from the whole discussion:

Does any one see's what this is happening and why it holds the equation ?
Or anything else I'm missing ?


texts: What's wrong with my M-ary ASK modulation demodulation MATLAB code?

texts:
My MATLAB code linked below implements ASK modulation of an M-ary signal, impresses it on a carrier, adds AWG noise of various SNRs. Then, for each SNR, demodulates the signal and calculates the Symbol Error Rate (SER). Finally, it plots  two $10\log_{10}(E_{bav}/N_0)$ vs. SER graphs, as shown in the image. One of the graphs (red) is obtained from the the theoretically calculated SER from the equation,
$$\text{SER} = \frac{2(M-1)}{M}Q\left( \sqrt{\frac{6 \log_2(M)E_{bav}}{(M^2-1)N_0}} \right)$$ $E_{bav}$ being the average energy per bit and $N_0/2$ being the variance of the added noise. Another graph (blue) is obtained by plotting the SER obtained from the code's demodulation scheme.

My questions are

What's wrong with my code that the two graphs are way apart?
Why even the theoretical curve does not give accurate values? (I checked it from a book)
Is my calculation of $N_0$ (Line 41, 42) correct?
Is my calculation of the basis function $\psi(t)$ (Line 35) and the expected cross correlator output (Line 36) correct? I calculated $E_g=T_s$ in these cases from
$$ E_g = \int_0^{T_s}g_T^2(t) dt = \int_0^{T_s}dt=T_s$$ by assuming the arbitrary pulse $g_T(t)$ based on which the the modulated waveforms are created is given by
$$g_T(t)= \begin{cases} 1 ~\text{for}~ 0 \le t \le T_s\\
 0 ~\text{otherwise}\end{cases}$$ Wrong assumption ?

Code link :
https://paste.ubuntu.com/p/jTs672h5z7/


texts:
I'm new to FFT and DSP and I want to ask you some questions about calculating FFT in Matlab. The following code is from Matlab help, I just removed the noise.

What happens if I choose the length of signal L > NFFT? and what's about choosing L different form NFFT?
In the following code, I'm not sure if I used window correctly. But when I use window (hanning in the following code), I can't get the exact values of amplitudes?
When L and NFFT get different values (I mean when I change L and NFFT,) then the values of amplitudes were different too. How can I get the exact value of amplitude of input signal?

I thank you very much and look forward to hearing from you :)

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 512;                     % Length of signal
NFFT=1024;                   % number of fft points
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);    input signal
X = fft(hann(L).*x', NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(X(1:NFFT/2+1)))     % Plot single-sided amplitude spectrum.


texts: Is there a closed form solution for a temporal calibration of two different signals?

texts:
Given two signals $f_1(t),f_2(t+dt)$ where $t$ is time and $dt$ is the time delay between the two signals,

is there any closed-form solution with respect to $dt$?
what are the efficient global solutions to this problem?

In my problem, the function is $f\in SO(3)$ but I can reduce it to $f \in R^3$ or just look at only one dimension.


texts: Ideal high pass filter for discrete signal

texts:
there. I currently get stuck on a question. I was asking to find an inverse discrete-time Fourier transform for the ideal high pass filter.
Here is the question
It is getting more confused after I read the second solution of this post
How to produce a high-pass filter from a low-pass one?
It tends out I can either calculate it by time-shifting property

Or I can calculate it through the definition

But they lead to a different result. My textbook solution said the second one is correct. But I cannot find any problem with the first one.
If we assume $\ 0 \le n < N-1$ for $h[n]$, using N = 61, $\alpha =30$, cutoff analog frequency $=1500Hz$.
(Assume the sampling rate $=44100 Hz$ )
I got these two graphs:
(for first solution)
and

(for the second)
Here are my codes

I really don't know why they lead to different results.


texts: Stationarity, discrete-translation operator, and the power spectral density matrix

texts:
Let $\mathbf{T}$ be the translation operator/matrix in discrete-time domain which can be written as $\mathbf{T} = \mathbf{\Phi} \mathbf{P} \mathbf{\Phi}^*$ where $\mathbf{P} = \exp(-i Diag([w_0, w_1, \ldots, w_M]))$ is a diagonal matrix where $w_k= 2\pi(k-1)/M$. We know that the stochastic signal $\mathbf{x}$ is called wide-sense stationary (WSS) w.r.t. the translation operator on in discrete-time domain $\mathbf{T}$ if and only if for all $r$:
${\mathbf R}_{\mathbf{x}} = \mathbb{E}[\mathbf{x}\mathbf{x}^{*}] = \mathbb{E}\big[(\mathbf{T}^{r} \mathbf{x})\big(\mathbf{T}^{r} \mathbf{x}\big)^{*}\big].$
We need to show that: a process in second-order moment is WSS if and only if
\begin{align}
{\mathbf R}_{{\mathbf{x}}} = \mathbf{\Phi} {\mathbf S}_{{\mathbf{x}}} \mathbf{\Phi}^{*},
\end{align}
where $\mathbf{S}_{\mathbf x}$ is a diagonal matrix with non-negative entries on its main diagonal and $\mathbf{\Phi}^{*}$ is the DFT matrix. How can we prove this result?


texts:
I have a signal which is sampled in real-time. At every time interval I receive a value. I need to characterize this signal as active or inactive (to take priority in the control actions). By active it means that we observe big changes over a specific time-window.

One way to do it as I consider, is to keep a history of N samples and calculate the mean value and the standard deviation with an update algorithm. Based on the std I can say if it is inactive or not. The problem is that I have tens of thousands of signals and the sampling period is not the same on all of them. So, keeping a variable length sample history for all the signals will be difficult.

For the mean value I can use averaging with forgetting factor. The forgetting factor can be connected to the sampling period of each signal to have a similar observation window on all signals.

Now the question:

Is there any metric similar to averaging with forgetting factor that can give me the same information as std? It has to be low on memory usage and computationally light (since it's a real-time operation).

Any papers, books or references are welcome!

Thanks in advance.


texts: Difference between speech and music signal

texts:
What are the difference between speech and music signal? Can anyone help?

I am doing project in audio processing. In that project I want make decision, that is which type of signal is given as input whether it is speech signal or music signal. I searched a lot in internet. But i could not find the required details.

If anyone knows about this please share the difference between the speech and music signal.


texts: When to use Fourier, Laplace and Z transforms?

texts:
If we have an LTI system, with an input signal $x(t)$, impulse response $h(t)$ and output $y(t)$, I was under the assumption that if the input and impulse response were continuous in time, then you would use the FT on $x(t)$ and the Laplace transform on $h(t)$, and multiply them together to find the output $Y(f)$ in the frequency domain.

Likewise, if $x[n]$ was a discrete input signal, $h[n]$ was a discrete impulse response, and $y[n]$ a discrete output, then I though we would use the DFT on $x[n]$, z transform on $h[n]$, then multiply them together to get $Y(k)$ which is a discrete output in the frequency domain.

However I've seen people saying we would use the Fourier transform on the impulse response to get the transfer functions here so it's kinda thrown me off a little bit. Could anyone clarify this?


texts:
Images show a discussion I picked up from a PhD thesis about a cyclostationary process and need help interpreting it.

"In the time domain the upsampling process creates a signal whose distribution of amplitude is periodic or cyclic in nature" : I see what the upsampling has done but I fail to understand how the amplitude distribution is periodic?

The author goes on to say "yet within each of the  four  time  offsets  the  distribution  is  stationary,  leading  to  the  name cyclostationary" : is it stationary within the time offsets?

" In  the frequency domain four perfect images of the original spectrum are created due to the upsampling" how has upsampling caused this to happen?

I don't mean to discredit the authors work, this is merely a reflection of my lack of understanding.


texts: How to calculate x and y of any object inside an image when stretched?

texts:
Suppose a image has

Height = 200px;
Width = 100px;

A point in that image which has

Y axis = 40px
X axis = 30px

Now if I resize that image:

New Height = 300px;
New Width = 150px;

  Is there any formula to calculate new Y and X axis of that point?


texts: Acoustic waves propagation distance relation with carrier frequency?

texts:
In electromagnetic waves, the energy formula is $E=hf$, that means frequency increasing will increase energy and so on the propagation distance.
In underwater acoustic waves, low frequencies can travel larger distances comparing to higher frequencies.

The question is: isn't the relation $E=hf$ is applicable in underwater acoustics?
If it is applicable, does it mean that an increase of carrier frequency will increase the transmission distance of underwater acoustic waves?


texts: 8*8 block matrix in JPEG image compression?

texts:
In standard JPEG format of a image, Discrete Cosine transform is used. But instead of applying the transform on whole image, we first divide the image in 8x8 block and apply transform on each of them. Thus during quantization we remove small coefficient in higher frequency. These steps are explained in details here.

But as we know that in N-samples DCT we can have only N frequecy coefficient. If we apply the DCT on the whole image we will able to get much higher frequency coefficient.

Can someone here highlight any disadvantage of applying DCT over the whole image at a time?


texts:
My question is, what frequency ranges do MFCCs correspond to?
My line of thinking thus far is that if your highest frequency in the signal was 8000 Hz, you could convert to Mels (= 2834 Mels), take 13 (for 13 MFCCs) evenly spaced measurements (0, 118.12, 236.24, etc.) and then convert back to Hertz to get start and end points for where the triangular windowing functions would be. Does that mean that MFCC1 corresponds to 0 to 77.5 Hz, MFCC2 corresponds to 77.5 to 163.5 Hz (approx, I know there is some rounding here)? How could I visualize these triangular windows to see what frequency value had the most weight (since I guess the frequencies further from the middle of the triangle would be represented less)?


texts:
I want to design a FIR filter by using  'cfirpm' function , which can i use the group delay and grid frequency,
my code is :

N=30 ;
f=[linspace(-1,-.5,120),linspace(-0.4,0.7,60),linspace(0.8,1,90)];
d=[zeros(1,120),ones(1,60),zeros(1,90)].*exp(-j*pi*f*N-1/2);
b=cfirpm(N,d,@lowpass);
fvtool(b,1,'OverlayedAnalysis','phase');

But i got the following error in matlab,

Error using cfirpm
Band edges must be monotonically increasing.

How can i solve this problem? or,how can i add the group delay in this function?

Thank you in advance.


texts:
As a newbie I am trying to understand I and Q modulation and demodulation. So from what I read the I and Q components of a signal separately mix with RF signals that are 90 degrees out of phase with each other and then the summation is transmitted as an RF signal.

So let's say I am listening to a voice of a reporter via FM radio. The reporter's sound wave gets converted to electricity with a microphone and then presumably gets converted to EM wave with an antenna. I am kind of lost where modulation happens and how would his sound voice gets divided to I and Q before modulation. Is it done mathematically and digitally?


texts: FSK carrier frequency estimation in non equiprobable data

texts:
We know that in equiprobable data,change in phase over a time interval, this result will be directly proportional to frequency:
f∝Im(s(n)s(n−1)∗)
Where s(n) is the normalized complex FSK time domain signal.

In the MPSK Modulations, Cyclic Spectrum Density method can give a peaks, for example for the QPSK:

But in the in the non equiprobable FSK signal, this peaks do'nt appear. How to solve it?


texts: High pass filter of 0.01 Hz for EEG slow drifts

texts:
I am trying to make a high pass filter for electroencephalographic (EEG) datasets, in order to get rid of very slow drifts. However frequencies around 0.3 Hz are very important for the research these data serve.

I use Matlab and I tried various things. It seems that a very low dB attenuation, of 3 dB is fine with FIR filters. For example I used that:

h=fdesign.highpass('Fst,Fp,Ast,Ap',0.005,0.01,3,1,250); % i needed to cut
d=design(h,'butter');
fvtool(d)

However this filter makes a dc drift at the start of the data and I would not wish to loose these first datapoints.

I have read through your pages that a clever idea is to design a low pass filter that later I would subtract from the real data. I did it with an equiripple filter and it worked but left a baseline lift. I do not want to run now average removal as this is in a next step in the analysis protocol.

Any suggestions?


texts: Maximum Likelihood Estimator (MLE), MMSE and LS - Are All of Them Regressor, Estimator and Predictor?

texts:
Can all three criteria ML, MMSE, and LS be called regressor, estimator, and predictor ?

If not, an intuitive explanation of why they can't be, would be good.


texts: changing circular shift in hologram reconstruction

texts:
I'm trying to write a very simple code for the reconstruction of an in-line hologram of particles in a micro channel (digital holographic microscopy), but I'm getting some results that are hard for me to understand.

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

hologram = imread([file_path file_name]);
N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win');

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram
fft_shifted_hologram = fftshift(fft_hologram); %Fourier shift of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram
psi = G .* fft_shifted_hologram;

% get slice of reconstructed optical field through inverse FT
recon_image = ifftshift(ifft2(psi));

%write image to file
imwrite(uint8(abs(recon_image)), [file_path_save 'out' num2str(k,'%2.2d') '.tif']);

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(-N/2:N/2-1,-M/2:M/2-1);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
    ./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
    ./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root<0) = 0;

G = exp(-2*pi*1i*z/lambda * sqrt(arg_root));

The actual reconstruction appears to be correct, but I am getting a circular shift that varies with the reconstruction (propagation distance).
Here's the original hologram:

Here's z = 0.01, which (just) looks like I missed the ifftshift():

Here's at z = 0.0105:

At z = 0.011:

The circular shift keeps getting smaller, and by the last propagation distance (z = 0.2) it is practically gone:

I've looked at the code for quite some time and can't figure out what's causing the circular shift, so I would greatly appreciate it if anyone could explain to me what is causing/how to correct it. Thanks.


texts: Bessel filter second-order sections Q and Fc multiplier derivation

texts:
A practical guide for digital IIR audio filters has cookbook-style values for creating higher-order Bessel filters out of biquads, but the values listed aren't very precise:

You should multiply the Fc for each stage by the following coefficients:
1: 1.00
2: 1.27
3: 1.32 1.45
4: 1.60 1.43
5: 1.50 1.76 1.56
6: 1.90 1.69 1.60

The corresponding Q values are:
1: ---
2: 0.58
3: --- 0.69
4: 0.81 0.52
5: ---- 0.92 0.56
6: 1.02 0.61 0.51

What are the expressions to calculate these values exactly?

For the Butterworth table, the values are given by Q = 1/(2 * sin((pi / N) * (n + 1 / 2))), for instance.  So where the table says 4: 0.54 1.31, the equation gives 1.3065629648763766 and 0.54119610014619701.  But I can't find the expression for Bessel filters, except that 2: 0.58 is $1\over\sqrt3$.

Prototype-generating functions just have long tables of numbers, so maybe these are not easy to calculate?  (Bessel functions?)  If so, a listing of these Q and Fc values for the first several orders would be good.

(And yes, I know the bilinear transformed Bessel is no longer linear-phase.  :/  I just want to know where these numbers are from.)


texts: What Is the Difference Between the Terms Color Intensity and Color Saturation?

texts:
In Digital Image Processing, what is the difference between the intensity and the saturation of a color?

How would you explain the difference to a newbie?


texts: I am working on tracking an object in a video by using its color. what i did is that i have converted the color of a selected object into LAB color space since it is a uniform color space. and then in each frame of video i am searching for same color by finding color difference in LAB color space. Now main problem is arising when more than one object of same color is appearing in the video scene. one method i adopted for avoiding detection of wrong object is to search for same color in limited area  around last detected position of the object. but this is causing failure to detect fast moving objects and leaving my methods suitable for slow moving targets only. so i want to know is there any method i can apply for avoiding wrong detection.


texts: Bose Chaudhuri Hocquenghem (BCH) vs Reed Solomon (RS) Error correcting codes and sound transmission channel modeling

texts:
I've implemented an acoustic based information transmission system using a bandwidth of 2 kHz and the band between 18.5 and 20.5 kHz. I'm trying to find the best error correcting code for it, in the sense that I could get the same error correction capability with the lowest redundancy added (I want to get the best speed possible). I have been able to transmit data between 2 smartphones at low rate (400 bps) and short distances (5-10 cm). At the moment I've design a package that I send between peers and I'm using 4 BCH(63,39,4) blocks there. I have also added some parity check bits (just for error detection) in case BCH fixed the incoming package in the wrong way (could happen if more that 4 errors occur in any of my 4 BCH blocks). I interleave the BCH coded bits before I modulate them trying to equally distribute the error burst among the 4 BCH blocks. Right now I'm capable of implementing both RS and BCH coding schemes, but I was wondering if there was some criteria I could follow in advance for choosing which one might perform better, would it be worth the time implementing RS. I read this document which says that for a Rayleigh Fading Channel being simulated in Matlab, BCH outperformed binary RS for a specific N,K,T selection. But they mention that this has a strong channel model dependency. My questions are:

What channel model fits data transmission over air best, I guess distortion generated by speakers should also be take under account, specially at the frequencies at which I'm working?
Is there a general study comparing BCH vs RS in different channels?
Is RS theoretically better than BCH for the same code rate n/k. Non binary RS may be considered, I'm just interested in same error correcting capability with greatest bit-rate ?

Independent answers to any of these questions are most welcome

Note: At this low bitrate CPU is not a bottleneck. Anyway the most time consuming task of my system is FFT based matched filter demodulation (specially if I use high order modulations), not BCH error decoding with such small code length.


texts: Low pass butterworth filter is distorting

texts:
my low pass filter keeps distorting and I'm at a loss to explain why. It calculates the correct coefficients, but distorts everything. The code is as follows. It's driving me mad. What's weird is my high pass works perfectly and it's coded in a similar way, but with a negative b1 coefficent and the same a1 and a2 coeffcients.

float cutoff = 100.00f * M_PI * 2.0f;
float c      = 2.0f * sampleRate;
float q      = sqrt(2.0f)/2.0f;

g_a1LP = (2.0f * pow(c, 2.0f) - (2.0f * pow(cutoff, 2.0f))) / (pow(c, 2.0f) + ((cutoff*c)/q) + pow(cutoff, 2.0f));
g_a2LP = -(pow(c, 2.0f)  + pow(cutoff, 2.0f) - (cutoff*c)/q ) / (pow(c, 2.0f) + ((cutoff*c)/q) + pow(cutoff, 2.0f));
g_b0LP = 1.0f;
g_b1LP = 2.0f;
g_b2LP = 1.0f;

The values for a1 = 1.97985 and a2 = -0.98005. The actual processing is then

//Low pass
            g_lastLowPassOutSamples[0] = (g_b2LP * g_lastLowPassInSamples[2]) + (g_b1LP * g_lastLowPassInSamples[1]) + (g_b0LP * g_lastLowPassInSamples[0]) +  (g_a2LP * g_lastLowPassOutSamples[2]) + (g_a1LP * g_lastLowPassOutSamples[1]);


texts:
I know all about conversion of frequency from HZ to normalized frequency and from normalized frequency to angular frequency ,but I wonder how it is done in MATLAB.
I have the following code:

f:0:0.001:1;
omega=pi*f

here f=normalized frequency and omega=angular frequency. My question is how do we get the interval value(0.001) for our normalized frequncy.

Thanks for any help


texts: Confusion Related to the Laplacian Filter

texts:
I have this confusion related to Laplacian filter. It states that it uses second derivative and finds the point of zero crossing(where the second derivative is zero) for which the first derivative has the high peak. I have seen them using kernels like

0 1 0
1 -4 1
0 1  0

I am not sure how this kernel finds the point of zero crossing. I mean definitely this kernel is equivalent to the second derivative.

Suppose I have points something like this

0 255 255
0 255 255
0 255 255

Definitely, it has an edge at the middle. But if I convolve the kernel with this image for example. I won't get zero result. So how come its called zero crossing


texts:
is there a way to just perform block motion compensation using ffmpeg command line (i.e. no addition of residual encodings). I want to compare the quality of the video frames produced using various block-motion algorithms.

I want to be able to extract motion vectors and generate motion compensated frames as in the SckitVideo Library (http://www.scikit-video.org/stable/motion.htmlhttp://www.scikit-video.org/stable/motion.html). This code is however limited to I,P,P frame sequences.
There must be a way of achieving the same result for I, B, B, P, B, B, I sequences.

The extract_mvs.c example in FFmpeg's github repository uses FFmpeg's libav library to extract motion vectors from side information in the video encoding.(https://github.com/FFmpeg/FFmpeg/tree/master/doc/examples)

At the moment I am looking for a way to do the compensation (ffmpeg command line would be preferable but reference to coding examples using ffmpeg's libav library would also be much appreciated).

Thanks in advance!


texts: Demodulation of OFDM symbol in LTE for multi-user case

texts:
I was trying to understand what exactly happens when you try to modulate an ofdm symbol in the downlink in LTE for multi-user case.

When we modulate, usually we allocate different sets of RBs for different users. Since an RB is 12 subcarriers in frequency, a given ofdm symbol (say 2048 SCs in a 20 MHz case) can contain data from different users. And we do take IFFT of this 2048 sized data vector comprising of different user data. Is that right?. OR do we append 2048-12 zeros to a specific user data and then take IFFT?

Thanks in advance.


texts: Peaking filter (loudness) cutoff frequency in an equalizer

texts:
An equalizer we have consists of Loudness ,Bass ,Mid and Treble .The loudness is nothing but a Peaking filter with cutoff frequency 20Hz .I read this in a book that outlines the EQ design.
I was wondering why the Loudness block which is basically a Peaking Filter should have a cut off frequency of 20Hz ?? Any intuition behind this


texts: MATLAB: Why my Pixel Indices are so different from that in Spatial Coordinates

texts:
The question is revised to make it more specific and constructive.

I am doing imtransform, yet the xdata and ydata are very much different from that of size(output).

[ output, xdata, ydata ] = imtransform( imgDouble{1}, tform );

Real data test:

K>> size(output)

ans =

         854        1230

K>> xdata(2) - xdata(1)

ans =

   1.6635e+03

K>> ydata(2) - ydata(1)

ans =

  876.1373


texts: Measuring signal phase difference between receivers

texts:
I am working on a project which aims to locate a transmitter.

I consider doing this using several receivers and measuring the phase differnces between them.

My problem is that arduino measures once every ca $dt=10^{-6}$ seconds, and therefore, to notice a phase difference, the distance between two receivers must be at least $c\cdot dt$ which is ~300 m.

I know that there are methods of finding the phase difference in much smaller systems, but I can't figure out how to solve that problem.

So, how can I measure the phase difference without placing the receivers that far away from one another?


texts: Determinant of Hessian approximation (SURF)

texts:
I have a question regarding formula in SURF article by Bay et al.

Theory

Given a point $p=(x,y)$ in an image $I$, the Hessian matrix $\mathcal{H}$ in $x$ at scale $\sigma$ is defined as follows

$$
\mathcal{H}(p, \sigma) = \begin{bmatrix}L_{xx}(p, \sigma) & L_{xy}(p, \sigma)\\L_{xy}(p, \sigma)& L_{yy}(p, \sigma)\end{bmatrix},
$$

where $L_{xx}(p,\sigma)$ is the convolution of the Gaussian second
order derivative $\frac{\partial^2}{\partial x^2}g(\sigma)$ with the image $I$ in point p, and similarly for $L_{xy}(p,\sigma), L_{yy}(p,\sigma)$.

The $9\times9$ box filters in Fig. 2 are approximations of a
Gaussian with $\sigma=1.2$ and represent the lowest scale (i.e.
highest spatial resolution) for computing the blob response
maps. We will denote them by $D_{xx}$, $D_{yy}$, and $D_{xy}$. The
weights applied to the rectangular regions are kept simple
for computational efficiency. This yields
$$
\det(\mathcal{H}_{approx}) = D_{xx}D_{yy} - (wD_{xy})^2.
$$

The relative weight $w$ of the filter responses is used to balance
the expression for the Hessian’s determinant. This is
needed for the energy conservation between the Gaussian
kernels and the approximated Gaussian kernels,
$$
w = \frac{\big\lvert L_{xy}(1.2)\big\rvert_F\big\lvert D_{yy}(9)\big\rvert_F}{\big\lvert L_{yy}(1.2)\big\rvert_F\big\lvert D_{xy}(9)\big\rvert_F} = 0.912\ldots \approx 0.9,
$$

where $|X|_F$ is the Frobenius norm. Notice that for theoretical
correctness, the weighting changes depending on the
scale. In practice, we keep this factor constant, as this did
not have a significant impact on the results in our
experiments.

Questions

Main question is how this formula for $w$ is obtained?
And my another concern is to why $w$ is within the brackets with $D_{xy}$. If it's just a balancing coefficient wouldn't it make more sense to write like this $wD_{xy}^2$? So I guess there is some motivation behind it.

My general purpose is to generalize this method for 3D, so it would be great if anyone could share some thoughts about balancing coefficients in that case or some useful links containing relevant information.


texts: How to measure quality of an audio/music recording compared to the original

texts:
I have 3 wav files:

The original recording of a viola note
A noisy re-recording of the original clip
A different noisy re-recording of the original clip

I'm looking for a way to show which of the two noisy re-recordings is better than the other, and/or more similar to the original.

I can plot the time-domain waveforms and spectrograms which show that the signals are vaguely similar, but ideally I could find a measurement that gives me a single number - something like "recording 1 = 35% quality, recording 2 = 63% quality". I've tried "corcoeff".

In the image below, the original is in blue, the decent recording is in green, and the bad one is in red.


texts:
(Schuam's DSP Outline, 2nd edition, problem 6.11(c), page 241).

Is there a DFT down-sampling property that looks something like this:

Given $x[\![n]\!]_M$ we want to downsample from M to N to obtain $y[\![n]\!]_N$...

$$
\Bigg(y[\![n]\!]_N = x\Big[\!\Big[\  [\![n]\!]_M\  \Big]\!\Big]_N \Bigg) \xrightarrow{DFT} \left( Y[\![k[\!]_N = X\left[\!\left[ \frac{M}{N}k  \right]\!\right]_M \right)
$$

$$
\text{for } (M > N)\ \&\ (M \ modulus\  N = 0)
$$

For Example, if I downsample $x[\![n]\!]_6$ which was modulo 6 and convert it to modulo 3 and assign it to $y[\![n]\!]_3$ which is modulo 3, then it looks like this in frequency domain:

$$
\Bigg(y[\![n]\!]_3 = x\Big[\!\Big[\ [\![n]\!]_6\ \Big]\!\Big]_3\Bigg) \xrightarrow{DFT} (Y[\![ k]\!]_3 = X[\![ 2k]\!]_6)
$$

Is this relationship true, and is there an easier way to write this relation as a transform?

Maybe I can invent a math notation to keep my DFT table from getting to crazy:

$$
\Bigg(y[\![n]\!]_N = x[\![n]\!]_{M\ to\ N} \Bigg) \xrightarrow{DFT} \left( Y[\![k[\!]_N = X\left[\!\left[ \frac{M}{N}k  \right]\!\right]_M \right)
$$


texts: How to filter a signal sampled at different frequencies

texts:
I have a signal which was sampled using different frequencies. The 3 regions are:

0-10 ms: 10000 Hz
10-50 ms: 1000 Hz
50-200 ms: 100 Hz

The signal is very noisy, so I thought about smoothing it with a low pass filter. I worked with each of the 3 regions independently:

Apply window in time domain (Hann window)
FFT
Apply half Hann window, to remove high frequencies
IFFT

Each of my regions looks better now but I am having problems stitching them back together, due to the beginning and end of each stage deviating considerably from the original values.

What is the standard procedure of applying a low pass filter to a signal sampled at different frequencies during several stages?

Edit:

In case it might be relevant, this is the half Hann window used as a low-pass filter. The window length (in the example, 50) is different for every region. I am not interested in a specific cut-off frequency; the filter is used to remove noise and I am evaluating the output visually.


texts: Easiest way to plot phase/amplitude/impulse response for zeros/poles in z-plane

texts:
I want to experiment to understand how poles/zeros affect the impulse/amplitude/phase response.
So I wonder if there is any software where you can mark poles/zeros and get this? (I have not found any yet)
Or otherwise how can this easily be done in matlab?
For example
How can I get the responses then? o is zero and x is a pole.


texts: Estimating the bandwidth of various signals

texts: I am trying to estimate the bandwidth of different signals, with different modulations, as seen in a PSD graph.  I've heard about a "$3\ \rm dB$ rule" in this context.  I don't need extreme accuracy.  Might drawing the edge of the channels at $3\ \rm dB$ down from the top work for me?


texts:
While studying Laplace Transform i found that region of convergence (ROC) is important because for some problems we have same Laplace Transform but different ROC helps us to take correct inverse Laplace Transform.
So, Now i am practicing inverse Laplace Transform problem i found almost every probelm to find $x(t)$ is without ROC.
So i want to ask that without have ROC how can we solve inverse Laplace Transform?


texts: How can we evaluate the matched filter's impulse response of a Gaussian function $x(t) = \exp(-\frac{t^2}{2})$. As far as we know, for a signal of finite duration $T$, the impulse response of it's matched filter is $h(t) = x(T-t)$. But as Gaussian function is of infinite duration, how can we find the matched filter ?


texts: FFTW audio artifacts when modifying magnitudes in frequency domain

texts:
I'm currently working with FFT (FFTW3) that I'm using to apply treatments on audio files frequencies.
The Forward/Backward test is passed, since I can get the exact  same soundfile when processing IFFT. So I assume my implementation is clean.
The thing for which I'm not an expert is to treat data then. For example, I have tried to implement a really naive filter to remove some of the frequencies in the sound. It is simply written like this :
int fft_buf_size = fft_len / 2 + 1
for(int i = 0; i < fft_buf_size; i++)
{
  int freq = i * (sample_rate / fft_len);
  if(freq > 50 && freq < 400) {
    fft_buf[i][0] = 0.0;
    fft_buf[i][1] = 0.0;
  }
}

This seems to process successfully the treatment, buf it does generate artifacts after IFFT.
I have tried windowing samples before FFT and unwindowing after FFT (which seems pretty useless) with Hamming method.
Last thing I tried is to overlap my windows (with hop factor of 2 or 3), which seemed to me to be the best way not to produce those artifacts. But, it only changes the artifacts periods, as those artifacts are correlated to the fft length (window size) too.
The overlap/windowing system is done like this :

First reading

 //Before FFT, read samples

 // Move read cursor by fft_len / hop_factor
 sf_seek(infile, passed, SF_SEEK_SET);
 sf_readf_double(infile, buffer, fft_len);
 passed += fft_len / hop_factor;

FFT/IFFT

// initialisation of the fftw plans
 fftw_plan forward = fftw_plan_dft_r2c_1d(fft_len, inbuf, offt, FFTW_ESTIMATE);
 fftw_plan backward = fftw_plan_dft_c2r_1d(fft_len, ifft, obuf, FFTW_ESTIMATE);

// After each read with sndfile, I process buffers by channels
for(int ch = 0; ch < number_of_channels; ch++) {
  for(int i = 0; i < fft_len; i++)
     inbuf[i] = buffer[i * number_of_channels + ch] * hamming(i, fft_len);
  fftw_execute(forward);
  // process my offt buffer with filtering or something else
  // Then copy it in ifft buffer

  // IFFT
  fftw_execute(backward);
  for(int i = 0; i < fft_len; i++)
   buffer[i * number_of_channels + ch] = outbuf[i] / hamming(i, fft_len);
}

Then writing

 //After IFFT, write samples
 int overlap_size = (fft_len * number_of_channels) / hop_factor;
 int overlap_index = (fft_len * number_of_channels) - overlap_size;

  // sum with overlap buffer
 for(int i = 0; i < overlap_size; i++) {
  buffer[i] += overlap_buffer[overlap_index + i];
 }
 // then copy buffer for next overlapping
 std::copy(buffer.begin(), buffer.end(), overlap_buffer.data());

  // Move write cursor to read_cursor position
 sf_seek(outfile, passed, SF_SEEK_SET);
 sf_write_double(outfile, buffer, fft_len * number_of_channels);

Am I doing something wrong, or didn't I figure the way to avoid those artefacts ?
I'm still not an FFT expert, so some of the common uses may be missing.
Thanks !


texts: How to show Image FFT with origin as DC?

texts:
The page Here explains how to do the discrete fourier transform on a 2d image.  It mentions this:

  In most implementations the Fourier image is shifted in such a way
  that the DC-value (i.e. the image mean) F(0,0) is displayed in the
  center of the image. The further away from the center an image point
  is, the higher is its corresponding frequency.

I have the implementation working, but I'm not sure how to shift my image to make DC be the center of the image.  Does anyone know how to do that?

The best I can think of is to push everything right 1/2 of width and up 1/2 of height, and wrap the pixels around, but this doesn't seem like the correct thing to do.

Does anyone know how to make this transform?  Right now, DC is in the lower left of my image, at pixel (0,0)


texts: Fourier Transform negative amplitude meaning

texts:
I am reading this example http://www.thefouriertransform.com/pairs/truncatedCosine.php

What does it mean to have some of the frequency components be negative in its amplitude ? I am not talking about the negative frequencies.


texts: Physical Interpretation of $\mathbf A$ beamforming Matrix

texts:
The directivity of an array is defined as:

$$ D = \frac{1}{2} \int_{-1}^{1} \lvert B_{u}(u)\rvert^{2} du$$

Where $B_{u}(\cdot) = \mathbf{w}^{H} \mathbf{v}(\cdot) $ is the beampattern of the array in u-space, $\mathbf{w}$ are the array weights and $\mathbf{v}$ is the array response vector. Expanding this expression, we get:

\begin{align} D &= \frac{1}{2} \int_{-1}^{1} \lvert\mathbf{w}^{H} \mathbf{v}(u)\rvert^{2} du\\
&= \frac{1}{2} \int_{-1}^{1} \mathbf{w}^{H} \mathbf{v}(u) \mathbf{v}^{H}(u) \mathbf{w} du\\
&= \frac{1}{2} \mathbf{w}^{H} \Big( \int_{-1}^{1}  \mathbf{v}(u) \mathbf{v}^{H}(u) du \Big) \mathbf{w}
\end{align}
Now we let:

$$ \mathbf{A} = \frac{1}{2} \int_{-1}^{1} \mathbf{v}(u) \mathbf{v}^{H}(u) du $$

How can one interpret the physical significance of the matrix $\mathbf A$?
Is it sort of an 'array response matrix', or something to do with spatial sampling?

NOTE: $\mathbf A_{mn} = {\rm sinc}\left( \frac{2 \pi}{\lambda}\lvert p_{m} - p_{n}\rvert \right)$, Where $\lambda$ is the wavelength of the incoming signal and $p_{i}$ is the position of the $i$th sensor.


texts:
This is the linear phase condition for FIR filters as expressed by my prof:

I don't understand why $G(f)$ can be negative. Isn't the Fourier transform expressed in polar form ? So the magnitude is always positive.


texts: When i want to see heaviside or unit step function in matlab, as shown in attached figure. i get first value of 0.5,why 0.5? why not 1,as unit step is thought to be ?


texts: What is Finite Rate of Innovation Signal?

texts:
I have read about Finite Rate of Innovation signal by Martin Vetterli in here. But i do not understand several basic things

The paper said that finite rate of innovation is the number of degree of freedom of a signal. But i cannot imagine a degree of freedom in signal term. If i have a robot that move lateral in X, Y, and Z, hence i can say it has 3 degree of freedom. How can i imagine it in signal term?
A standard sampling is just sample a signal for every time interval, but in FRI sampling, there is a smoothing kernel, what is the purpose of this?
what is annihilation filter?

Thank you


texts: how to do initialization of random vectors in matlab

texts:
I have to generate random vectors for $\mathbf{p}^{(0)}$, $\mathbf{u}_l^{(0)}$,  and $\mathbf{v}_l^{(0)}$ for $l=1\ldots L$.
Now I am confused about how to generate these vectors as $l=1\ldots L$ is the number of users, for example, I have $L=2$ then the generation of vectors should like $\mathbf{p}^{(1)}$, $\mathbf{p}^{(2)}$ for $\mathbf{u}_1^{(0)}$, $\mathbf{u}_2^{(0)}$, $\mathbf{v}_1^{(0)}$, $\mathbf{v}_2^{(0)}$ and after that, I have performed steps 1 and 2 respectively, for each $l$. Kindly guide me I am very confused.
L=4;  N=4;  N_l=2;  P=10;

w = ones(L,1);
beta = ones(L,1);
%p = zeros(L,1);
H = sqrt(.5)*(randn(N_l,N)+1i*(randn(N_l,N)));
u_l = sqrt(.5)*(randn(N,1)+1i*(randn(N,1)));
u_i = sqrt(.5)*(randn(N,1)+1i*(randn(N,1)));
v_l = sqrt(.5)*(randn(N_l,1)+1i*(randn(N_l,1)));
u_l = u_l./norm(u_l,2);
v_l = v_l./norm(v_l,2);
for k = 1:L
    p(k,:)=P/L;
    w'*p(k)<=P;

A general code for initialized vectors are given in code but they are for $L=1$, I have work for generalized case where $L$ can be changed accordingly


texts: Why is there only one integration in the solution if there is two integral in the formula?

texts: In this problem the random variable is theta and according to the formula there should be two integrations but in the solution there is only one . Nor am i able to understand the meaning of x1 and x2 in the formula. this problem is driving me nuts..


texts: What is the best way(In terms of speed) to implement a direct form type 2 IIR Filter in VHDL?

texts:
I am trying to implement a 4th order direct form type II IIR filter in VHDL(Specifically, in xilinx FPGA) and when i implement it i get correct output yet very low maximum clock rate(Around 80 Mhz).
So i wanted to know what is the best way to implement IIR filters in terms of speed?Or at least where can i read more about IIR filter implementations in VHDL?


texts: Recover myself from Dilate and Erode for noise reduction

texts:
i have obtained some skin detection map (i post below two example of detection map). but, i'm using a normal webcam and the map is noisy. Not a lot, just a bit.
I think i can remove noise only with some combination of dilate and erode, and then threshold it.
But, like a stupid child, i have spent hours in trying to set parameters, use in different combination like: dilate twice, erode three times, blur, and an erode again and so on.
The only result i have obtained is to lost myself: it is not a robust approach at all.
So, for image like this, how can i use those super-popular filters? Or some advice for getting a nice contour of the face and removing the noise?

EDIT:
Those two are two different approach to skin recognition: backprojection and mahalonobis distance, using as training data the output of a simple haar face detection. I'm using both of them because with some illumination setting the first is better, for others the second is better. For my porpouse i can choose the best anytime I'm using my application.


texts: Resampling scattered data with MATLAB's $\tt interp$ and $\tt resample$

texts:
I am recording acceleration data with an MPU6050 connected to a Arduino1 and stored on a SD. Here you can find the code.

I need to calculate the FFT of an acceleration signal that was not sampled uniformly, so i have to resample my signal. After this suggestion I tried to study how to resample the signal.

From Mathworks documentation here it seems that I have to know the frequency of my signal and a nominal frequency. I read and studied the function resample, but I don’t understand how it works very well. In my specific case how to define

x= sin(2*pi*f*T)

But in my code below I tried to implement the function (to give you more information for a detailed answer and prove my efforts, but I don't have enough reputation to post all the graph).

Then I studied gettdata  “Scattered interpolant class” finding not very useful in my situation, they seems to work for 2 or more dimensions.

So I followed another way, I used the function interp1. I defined my equispaced vector using the mode and median value of my sample frequency, I discarded the average to avoid the influence of any peaks.

Is a correct approximation?
What kind of errors can occur?

This is my first time I approach to the FFT, and I have a civil engineer background, so is my first time with signal analysis, never studied Signal Theory before.

My interest for the FFT is to define the best low pass filter (example apply a Butterworth but I don’t know how to choose the filter order and cutoff frequency.

I don't have enough reputation to post more than 2 links sorry for that.

filename= uigetfile ('.txt');
fileID = fopen (filename);
logmpu6050 =csvread(filename);
fclose (fileID);
%Starting creating the specific Vectors
%Record Time in millisecond
time=logmpu6050(:,1);
%The x y z are converted to m/s^2
confactor=19.6/32768;
ax=logmpu6050(:,2);
confactor=19.6/32768;
ax=ax*confactor;
ay=logmpu6050(:,3);
ay=ay*confactor;
az=logmpu6050(:,4);
az=az*confactor;
%Define the sample rate subtracting from Sampletime i+1 Sampletime i
n=length(time);
for i = 1:n-1
samplerate(i) = time(i+1)-time(i);
end;

%I try to define the best frequency for "resample" function based on mode and median value
%I will not use average to avoid conditioning due to extreme value
f1=mode(samplerate);
f2=median(samplerate);
x1= sin(2*pi*f1*time);
x2= sin(2*pi*f2*time);
plot(time,x1,'.');
figure
plot(time,x2,'.');

%resampling at time mode record frequency
[test1,test2]= resample(az,time,f1);
%resampling at time median record frequency
[test3,test4]= resample(az,time,f2);
figure
plot(test2, test1)
figure
plot(test4,test3)

f1=int32(mode(samplerate));
f2=int32(median(samplerate));

%Try to implement the matlab function "interp1"
%My input are
%x-> time when the value is recorded "time"
%v-> the acceleration value "az"
%xq-> I tried two define two type of equispaced vectors
%xq1-> Defined by the total recording time divided by the mode
%xq2-> Defined by the total recording time divided by the median
totrectime= int32(logmpu6050(n,1)-logmpu6050(1,1)) %milliseconds
numsamples1= idivide(totrectime,f1); %Rounded toward zero to the nearest integers.
numsamples2= idivide(totrectime,f2);
f1double=double(f1);
f2double=double(f2);
for i = 1: numsamples1
xq1(i)= f1double*i;
end;
for i = 1:numsamples2
xq2(i)= f2double*i;
end;
xq1double=double(xq1);
xq2double=double(xq2);

vq1 = interp1(time,az,xq1double);
vq2 = interp1(time,az,xq2double);
figure
plot(xq1double,vq1);
figure
plot(xq2double,vq2);


texts: Why is a first/second derivative useful in spectroscopy?

texts:
I'm currently working on Raman spectroscopy and while reading some literature I came across the first and second derivative of a spectrum. It's not clear to me why they are useful to look at since they look very noisy to me.

I tried graphing the first and second derivative of a spectrum I attained from liquid whey (see the following pictures, 1. accquired raman spectrum (smoothed), 2. first derivative, 3. second derivative):

What information can I read from the first/second derivative in this specific example and in general?


texts:
I have big database of audio files, about 1-2 s long. Different people say different words, but during process I done some duplicates which I have to find and remove (database is about 100000 sounds and it is difficult to listen, have some noise, or are similar to human, but are different). Some duplicate audio files cant be shifted by half a second or less during cutting process.

Please help.
How can I find real duplicates in my audio database?


texts: Prevent negative image with BPSK modulation in a complex signal

texts:
In a baseband complex signal, I need to generate a BSPK modulated signal at a given frequency.  If I generate a CW at this frequency, all is good.  However, when I then BPSK modulate the CW, I get a signal at both the positive and negative frequencies.

Is there a way of implementing BPSK modulation without this negative frequency image, or must I simply implement a filter to remove the image?


texts: Is there a term called "Fourier Beamforming"?

texts:
I understand delay and sum beamforming, where we want to add the signals obtained from different sensors. But adding different time samples from different element of the array depending on the angle we are interested. For broadside or zero degree we would add the samples from the same time instant from each element of the array. On the other hand if we are interested to beamform in a particular angle there would be time delay between the samples of each element.

The question where "Fourier beamforming" comes into this picture? I have seen this term in many places. Any comment would be appreciated.


texts:
I have the following signal:

I'm trying to compute a Spectrogram algorithm, but, don't think I'm doing it right..

I have computed the following:

1) STFT (size 256 with an overlap of 128)
 2) Computed the logs using: '10 * log10(sqrt(re * re + im + im)

This is the result that I get:

But when I use pylab in Python (for the same signal):
  x = pl.specgram(signal)

I get the following result:

Using the matplotlib I get the following:

Obviously, these are very different results.. I don't know why I'm getting these, I'm new to signal processing and spectrograms. Hope someone can help

EDIT:

This is the result I have when doing imshow in python:


texts: Basic questions related to Kalman filter

texts:
Since I'm a newbie to Kalman filter, so, I had some confusions that I needed to clarify.

What is the difference between measurement and state?
Can they be of different dimensions?

All the examples that I have seen online
     take velocity and position into account while acceleration is
     external control variable. Now the state matrix is multiplied with
      a matrix containing $\Delta t$.

What happens to state matrix and the matrix it is multiplied with to get the updated state if my
  readings are acceleration itself? I am specifically referring to this dataset : Accelerometer dataset
Since readings already have
  noise, would we still add error/noise while calculating the
  updated state?


texts:
I have a signal in MATLAB

s = randn(1,1e6);

that is normalized in the following way:

s = s/sqrt(sum(abs(s .^2)) / 1e6);

I would be thankful if anyone could explain me why is the signal normalized in this way!


texts: How to Approximate / Fit a Function and a Lower Bound on It?

texts:
I am dealing with a complicated optimization problem with two parameters $N, K$. Setting $K =$ constant and finding the minimum for different $N$'s lead to a set of function values. I can do this for $K = 2,3,\cdots,8$ and $N = 5,6,\cdots,40$. Larger $K$'s or $N$'s take more time to calculate (e.g. 1 hour).

I've plotted what I get (blue) and a lower bound (green) for it, which is described by $L(N,K) = \sqrt{\frac{N-K}{K(N-1)}}$

$K = 3$

$K = 4$

Any of these two are desired :

An expression with respect to $(N,K)$ which best approximates the values I'm getting.
A tighter lower bound (for larger $N$'s)

What methods do you suggest?


texts:
I have been working with some seismic data and recently submitted a paper to a journal.

My seismic data has 493 traces: each trace are 8300 seconds long with sampling frequency of 2000 Hz.

One of the reviewers is asking me to generate a frequency spectrogram at the time= 2000s for all signals. My knowledge about the spectrogram is that you can only do it for one signal and there should be a time window to run the frequency analysis on. How can I do it?

Here are more detail about my problem:
These 493 traces recorded vibration around a cable in a oil producing horizontal well. The traces are 16ft away from each other along the well. The idea of the publication was to assess the spectrogram for traces  on a specific location that was the target for hydraulic fracturing. The figure below is the spectrogram of one of the traces (signals) close to the target zone. Note that there are 492 other traces and all of them have the same number of samples, frequency, and length. All these traces recorded vibration for the entire 8300secon. The question is how I could generate spectrogram of all 493 traces for example at 2000 second.
Here is the wording of the reviewer: "I welcome the addition of a spectrogram showing many channels at one particular time interval."
The paper is accepted, But I think it's better to address the reviewer's comment. I will appreciate your name in my acknowledgment.


texts: Reconstructing image from the dwt coefficients

texts:
I am trying to compress and reconstruct an image but facing problems with the data type of the compressed image pixel values which are of float64 and therefore, when I am trying to retrieve the pixel values from the compressed image , reconstruction is poor because on reading the pixel values corresponding to compressed image, the pixel value array elements are converted to uint8 type, thus leading to loss in the pixel data on automatic conversion from float to uint8. I want to take dwt of the image, store the dwt coefficients as an image (as images are of the lesser size as compared to pixel values stored in a file), and then retrieve the pixel values from the coefficient images, and thus reconstructing the image on taking inverse wavelet of the coefficients arrray obtained on reading the coefficient array.
I am implementing the complete system in OpenCV python.


texts: Bio-signal discontinuity after FIR filtering

texts:
I'm filtering respiration signal. I have a signal of duration equal to 10 minutes, so I decided to divided it into period of 10 s. I filtered each epoch with FIR filter, but I have discontinuities each 10s.
I attached an example.
I can not filter all the signal at once by the theory, so how can avoid discontinuity?


texts: Estimate Taylor series coefficients from samples of a function

texts:
Say I have measurements of a function $y = y(x)$, sampled at $x_i$ with some noise, that could be approximated by a Taylor series expansion. Is there an accepted way of estimating the coefficients for that expansion from my measurements?

I could fit the data to a polynomial, but that's not quite right, because for a Taylor series the approximation should get better the closer you are to a central point, say x = 0. Just fitting a polynomial treats every point equally.

I could also estimate the various orders of derivatives at my point of expansion, but then I need to make decisions about what differentiating filters to use and how many filter coefficients for each. Would the filters for the different derivatives need to fit together somehow?

So does anybody know of established methods for this? Explanations or references to papers would be appreciated.

CLARIFICATION

In response to the comment below, my sampling is a rectangular window from an infinite function, that isn't necessarily band-limited but does not have strong high-frequency components. To be more specific, I'm measuring the variance of an estimator (measuring displacement in a medical ultrasound signal) as a function of a parameter of the estimator (the level of deformation or strain of the underlying tissue). I have a theoretical Taylor series for the variance as a function of deformation, and would like to compare it to what I get from simulation.

A similar toy example might be: say you have a function like ln(x), sampled at intervals in x with some noise added. You don't know what function it really is and you want to estimate its Taylor series around x=5. So the function is smooth and slowly varying for a region around the point you're interested in (say 2 < x < 8), but is not necessarily nice outside of the region.

The answers have been helpful, and some kind of least-squares polynomial fit is probably the route to take. What would make an estimated Taylor series different from a normal polynomial fit, though, is that you should be able to shave off higher-order terms, and have the polynomial still approximate the original function, just within a smaller range about your initial point.

So maybe the approach would be to do a linear polynomial fit using only data close to the initial point, followed by a quadratic fit with a little more data, cubic using a little more than that, etc.


texts:
I'm trying to calculate the power of a signal and my tutor has given me this formula to do it

I've spent the past while building a program and now the foundations are there it's time to implement the maths side of it. The problem is I can't actually read it. Would someone be able to transcribe it and explain it to me?


texts:
I am working with speckle noise reduction in ultrasound images. I have used wavelet transform for removing the noise  by eliminating certain frequencies in order to eliminate any existing noise. Since in an image HH,LH,HL components contains most of the noise. We can eliminate noise by eliminating those components.

  step 1 : read an input image.
  step 2 : apply discrete wavelet transform.
  step 3 : Eliminate LH, HL, HH  by making these components zero.
  step 4 : Take inverse wavelet transform(let us consider this as output 1).
  step 5 : read the same input image and again take discrete wavelet transform.
  step 6 : eliminate LL component by making zero and keep LH,HL & HH because it contains edge information. (trying to remove noise from the edge information & add it back to Output 1)
  step 7 : take inverse wavelet transform .
  step 8 : Applying a Hybrid Filter (Median- Weiner) to the above image
  step 9 : We get the filtered image ( let us consider this as output 2)

I have represented the same in the below figure.

Now, I need to combine output 1 and output 2 . I just added both the outputs but am not getting proper output . Please suggest me the technique which I should use to combine these two output images. Kindly tell me whether it is possible or not and also tell me whether this concept is correct or wrong.


texts: Is interpolation of an audio signal to increase frequency resolution possible?

texts:
I apologize if some of what I ask is not entirely correct, I'm new to this field, but extremely interested.

I have an Audio signal of sample rate 44.1 kHz that I want to segment into 30 frames, and get the DFT of each frame to find the magnitude of certain frequencies for that frame. However, this means that I have a frequency resolution of 30 Hz in each bin, which isn't narrow enough.

Is it possible to interpolate the data to attain more data samples? As far as I'm aware, doubling the number of points would give a 88.2 kHz sampling rate, but still give a frequency resolution of 30. Would it be possible to treat the interpolated data as still having a sample rate of 44.1 kHz?


texts:
I am currently searching for methods of feature extraction from an ECG signal and I've stumbled upon the Karhunen–Loeve Transform. I've read some papers and I think I get the basics but my question though is, how do I extract the features(i.e R wave variability) based on the KLT transformed matrix?

Thank you in advance!


texts: How can I find SNR, PEAQ, and ODG values by comparing two audios?

texts:
I have two audios, A and B. B is the noisy version of A. I want to know how much B is different from A.
Can I do this by calculating signal to noise ratio (SNR), Objective difference grade (ODG), and Perceptual Evaluation of Audio Quality (PEAQ)?

How can I calculate them in Matlab?


texts: How to: pixel packs in an image to representative point coordinates (x,y)

texts:
Let's have an image (gray-scale or even binary) as shown on the following figure at the left hand side, the goal is to generate a list of points i.e., coordinates in the form of (x,y) for each pack of the dark pixels in the image.
What are the proper image processing tools to do this and where are they available?

Updates:
1)
Here you may find some more details about the problem. (Note the variation in the size of packs)

I may suggest having packs detected to compute the convex-hull boundary for each and then find the representative centroid {see this for details}.

2)
Here is the result produced by application of Distance Transform (suggested by "Libor"). Note my annotations on the figure. The method does not work as it was promising!

3)
Erosion eliminates small packs!

from __future__ import division
from scipy import zeros, ndimage as dsp
from pylab import subplot,plot,matshow,show

img = zeros((30,30))
img[10:14,10:14] = 1
img[16:17,16:17] = 1
img[19:23,19] = 1
img[19,19:23] = 1

subplot(221)
matshow(img,0)

subplot(222)
y = dsp.binary_erosion(img,[[1,1],[1,1]])
matshow(y,0)

subplot(223)
y = dsp.binary_erosion(img,[[0,1,0],[1,1,1],[0,1,0]])
matshow(y,0)

subplot(224)
y = dsp.binary_erosion(img,[[1,1,1],[1,1,1],[1,1,1]])
matshow(y,0)

show()

4)
Well here is a Python (i.e., the language of love :) ) implementation of the labeling idea (also proposed by "Jean-Yves" below):

subplot(221)
l,n = dsp.label(img)
sl = dsp.find_objects(l)
for s in sl:
    x = (s[1].start+s[1].stop-1)/2
    y = (s[0].start+s[0].stop-1)/2
    plot(x,y,'wo')

and the result:

Note that although it is done in Python so quick due to Scipy performance, the background procedure in label function should be an exhausting iteration. This may be considered as a trade-off. So for a while I keep being eager to seek more efficient algorithms. And also note that in the given code above I found the center of geometry so simply while for complex or asymmetric shapes this may cause positioning to be biased. That is it is a work in progress ;).

5)
Here is a complex case (a real image) captured from here on which the labelling proposal applied and you see the results. Note that it took only 0.015 s for whole procedure including labelling and finding the objects. Scipy guys, did very well job, I think. Wow! {right-click on image, click on view image for full resolution}


texts: Finding minimum-phase/allpass of transfer function

texts:
I am trying to find the minimum-phase system of the transfer function $H(z) = \frac{2 + 3.125z^{-2}}{1-0.9z^{-1}+0.81z^{-2}}$.

I know I need to find and remove the allpass (basically reflect the two zeroes on the imaginary axis inside the unit circle?), and I've gotten so far as this: $H(z) = 2\frac{1 + \frac{3.125}{2}z^{-2}}{1-0.9z^{-1}+0.81z^{-2}} = 2\frac{1 + \frac{3.125}{2}z^{-2}}{1-0.9z^{-1}+0.81z^{-2}} \frac{\frac{3.125}{2} + z^{-2}}{1 + \frac{3.125}{2}z^{-2}}=,$

and then I'm not really sure where to go next. Is my thinking correct thus far? What step should I do next?

Thank you.

Edit:
This is the kind of thing I'm trying to do, as per the example in my book:


texts: Signal compare - cross corelation in frequency domain matlab

texts:
I need compare two signals. Signals have simliar shapes but completely different values (as on picture):

So I think cross corelation in frequency domain should gave me a number that shows me how the signals are similar. But I have no idea how to do that in matlab.

Where I can find code to do it?


texts: On $H_\infty$ norm for transfer function

texts:
For a given scenario in the context of control system, I'm trying to investigate how the $H_\infty$ norm can be calculated for a transfer function as follows:
$$G(s)= \frac{w_n^2}{s^2 +2\zeta w_ns +w_n^2  }$$

$$\left \| G \right \| _\infty = \max \limits _{\omega} |G(j\omega)|$$

$$\left \| G \right \| _\infty = \sup \limits _{\omega} |G(j\omega)|$$

I know how to deal with such problems when the  case is specific one, for example the frequency value is known. However, I get confused with this general case. any idea


texts: Important differences between BCJR and SOVA algorithms?

texts: BCJR algorithm needs noise variance for Gamma (branch metrices) calculation. So I am wondering does SOVA (soft output Viterbi algorithm) needs noise variance?


texts:
Is it possible to write a FFT function, where there is fewer outputs than inputs? What in mean is, we would start computing 16-point FFT, and our result would end up be 12-point.
If it is possible than do we loose any data in the procoess of reduction?


texts:
I have recorded my own voice in Matlab and I intend to add some echo to it.I came up with one solution for getting the desired echo effect:

Delay the sampled audio in the time domain and adding it to the original sample.

In order to do this I'm using Matlab and I have basically done the following:

recObj = audiorecorder(44100, 16, 2);%sampling rate of 44100Hz, stereo
recordblocking(recObj,length);%record audio for a fixed length duration
y = time_delay(getaudiodata(recObj) , 5000 );%set a delay of 5000 to original sampled audio data

function [ y ]  = time_delay ( x , R )
%this function sets a delay to x of R
y = zeros(length(x) + R , 1);
y(R + 1:length(y)) = x(1:length(x));

As you can see I have a recObj wich I can use for playing the recorded audio.In order to get a delayed version of the recorded audio I have used getaudiodata(recObj) and my own function time_delay.A problem that I have is that I don't know how to get a new delayed audiorecorder object from the delayed audio data retrieved by the time_delay function and I need an audiorecord object cause this is the only way of playing something.
So, my question is:

Is this a good approach to get the desired echo effect in an audio
signal?
If it is, then how do I set the delayed audio data to an
audiorecorder object?


texts: Obtain a signal's peak value if it's frequency lies between two bin centers

texts:
Please suppose the following:

The frequency of a signal's fundamental has been estimated using FFT
and some frequency estimation methods and is lying between two bin
centers
The sampling frequency is fixed
Computational effort is not an issue

Knowing the frequency, what is the most accurate way to estimate the corresponding peak value of the signals fundamental?

One way might be to zero-pad the time signal to increase FFT resolution such that the bin center will be closer to the estimated frequency. In this scenario, one point I am not sure about is if I can zero-pad as much as I want or if there are some drawbacks in doing so. Another one is which bin center I should select after zero padding as the one I am obtaining the peak value from (because one may not hit the frequency of interest exactly, even after zero padding).

However, I am also wondering whether there is another method which may deliver better results, f.e. an estimator which uses the peak values of the surrounding two bin centers to estimate the peak value at the frequency of interest.


texts: Is it possible to recover the time-domain signal after manipulating it in the frequency-domain?

texts:
I've worked a fair amount with EEG signals, though I've never had formal training in signal processing, so please excuse my ignorance.

The problem is this: my signal has noise at many, many frequencies (6 Hz, 12 Hz, 38 Hz, ... 94 Hz, 100 Hz, 106 Hz...). See the following figure of the frequency-domain signal; the "real" signal is way at the bottom, not really visible at this scale. (Besides 100 Hz hum, I don't know what the sources of the noise are.)

% MATLAB code
L=2560;  %length of signal
Fs = 256;  %sampling frequency

Y = fft( signal);

% create one-sided frequency-amplitude figure
P2 = abs( Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2)) / L;

plot( f, P1)

What I would like to do is truncate all these noise peaks in the frequency domain, and then convert the signal back to the time domain, in effect filtering out all of those peaks. But I don't know if this makes sense, or is mathematically possible?

What I can do right now is, after the FFT, give each spectral component a random phase (following this post), thus "simulating" a signal with the given spectral properties.

% Assume I do the peak-truncation here
%%

% reverse from one-sided to two-sided
P1(2:end-1) = P1(2:end-1) / 2;
P2 = [P1(1:end-1) P1(end-1:-1:1)];

% add random phases
phi = unifrnd( 0, 2*pi, 1, length(P2));  %vector of random phases between 0 and 2*pi
Y_simulated = P2 .* exp( 1i * phi);  %the "simulated" FFT

signal_simulated = ifft( Y_simulated);

This produces a nice simulated signal which is similar to the original input signal, but is not time-accurate; you see that the peaks and troughs don't line up. (Note: This example is a clean signal and I did not do any truncation of noise peaks.)

I feel like there should be a way to take the phase information from the imaginary part of the FFT (Y) and reapply it at the end to recover a signal that is time-accurate. But how?

(Related question: Is this how filters work? Or do they work by a different principle?)


texts: Consider an ideal low pass filter $H(\omega)$, and the input to this filter is the periodic square wave $x(t)$. Find the output $y(t)$

texts: The solution to the problem is $$y(t) = 5 + \frac{20}{\pi} \sin(\pi t) + \frac{20}{3\pi} \sin(3 \pi t) $$ and to get that the solution says to find the Fourier series expansion of $x(t)$ and I am having trouble doing that. I find that $C_0 = 5$ and that $$C_k = \frac{5}{j k \pi} \left( 1 - e^{-j k \pi} \right)$$. I don't know how to compute $$C_0 + \sum_{k=1}^{\infty} C_k e^{j     k w t} $$. Also since the cut off is $4 \pi$ so when the Fourier series expansion greater $4 \pi$ those values will not be included as seen in the solution...


texts: Which Image pre-processing method is the best for Histopathology images? (H and E stained image)

texts:
I am working on classification of breast cancer histology images (Dataset:ICIAR BACH challenge 2018). But I don't really know which method is the best way to enhance the features of H&E stained image. I've tried the MATLAB code for Contrast limited adaptive histogram equalization method for enhancing the features, Kmeans Clustering, color normalization etc.
But nothing is giving a good result.
Can someone please suggest another method for preprocessing the H&E stained image for best result?
(FYI I am completely new in the field of image processing and machine learning. So I don't really have a clear idea in this field.)


texts: Feature Extraction of FFT for One Class SVM

texts:
I'm looking for a good way of extracting features from the frequency domain of vibration data for a one-class support vector machine.

The image below shows an example from the dataset found in this link.
The dominant peak found at around 960 Hz corresponds to the "main frequency" of the system. E.g rotation*number of gears.

Im keeping this question fairly open ended, because i know there are several methods for doing this and I'm mainly looking for your suggestions and experience.

Would it be possible to analyze this so I get some specific information regarding these peaks caused by the faulty gear tooth?


texts: FFmpeg audio filter pipeline for speech enhancement

texts:
Using ffmpeg, I'm recording a talk, its audio with a semi-distant microphone. It makes sense to somewhat enhance the speech clarity, reduce background static and environmental noise, etc., in the most generic terms.

What combination of ffmpeg audio filters am I best advised to use to enhance recorded speech?

I'm currently using

highpass-frequency=300,
lowpass-frequency=4000,

as that is the range speaking voice is most prominent in,

bass=frequency=100:gain=-50,
bandreject=frequency=200:width_type=h:width=200,

to additionally attenuate the bass frequencies, and

compand=attacks=.05:decays=.05:points=-90/-90 -70/-90 -15/-15 0/-10:soft-knee=6:volume=-70:gain=10

to apply dynamic range compression and supposedly make the softer passages (above -70dB) louder and the louder ones more even.

I understand digital filters are not magic, but I'm looking for something more from ffmpeg's (or LADSPA's) repertoire to enhance the sound more, make it more resilient to noise, sporadic coughs, perhaps lose the acquired metallic aftersound ...

The Linux command line I'm using to test:

filters='
...above filter lines, in sequence...
'

arecord -q -r 11025 -f s32_le -c 1 -t wav - |
    ffmpeg -f wav -i - -af "$filters" -f wav - |
    aplay


texts: Additive synthesis vs. "just draw the waveform"

texts:
I've been experimenting with PYO recently, a Python module for sound synthesis and audio processing.  This is really my first foray into this topic, so I have what is probably a very basic question.

There are a few alternative mechanisms available for creating periodic waveforms such as a sawtooth or triangle wave:

Using additive synthesis (via SawTable, which internally calls HarmTable)

This approximates a waveform through successive addition of sine waves.  E.g.,

from pyo import *
table = SawTable(order=20)

Which gets me:

Using a "linear" table (LinTable)

This allows one to construct a waveform as a series of straight line segments.  E.g.,

from pyo import *
table = pyo.LinTable([(0, -1), (8191, 1)])

Which gets me:

These two waveforms sound similar, of course, but they are qualitatively different.  What is the technical difference between the two?  Is there a reason to prefer one mechanism over the other?


texts:
there is a ODE
$dx/dt-ax(t)=u(t)$
Firstly, I have to find the transfer function $u$ to $x$ using exponential input. So I put $u(t)=e^{st}$ to the ODE, and I got $x(t)=x(0)e^{at}+1/(s-a)e^{st}$. In this case, can I say that $1/(s-a)$ is a transfer function?
Secondary, if $u(t)=e^{at}$, I have to show $x(t)=x(0)e^{at}+te^{at}$. I can derive x using undetermined coefficient method but how can I derive it using transfer function? Thank you.


texts: Goodness of fit for complex valued curves (i.e. frequency responses in frequency domain)

texts:
My apologies for perhaps the stupidity of this question.

Presume that one has the 'frequency response' $Y_{data}(k)$ of a system and also has an estimated model $Y_{syn}(k)$ that fits the data.

How can I quantify the 'goodness of fit' of $Y_{syn}(k)$ to $Y_{data}(k)$ ? Notice that both curves are complex-valued at every instant k (i.e. frequency domain).

For example, assume that I have two estimated models: $Y_{syn}^a(k)$ and $Y_{syn}^b(k)$. My goal (besides visualizing the residual) is to evaluate in a objective way the fit of both models to $Y_{data}(k)$ (the original data set). Is there a single metric to do this?


texts:
Usually in radio systems in order to move the frequency response of the message signal to a new frequency band centered at f, the message signal is mixed with $\cos(2\pi ft)$.

What would change in frequency if we multiply the message signal to  $\sin(2\pi ft)$ instead of  $\cos(2\pi ft)$?


texts:
I am trying to recover my baseband signal that was AM modulated. I am using the filter command in MATLAB to perform the low pass filtering ...the code is given below:

Fsig=10;                    Fc=100;              Fsamp=10.*Fc;
t = 0:1/Fsamp:2*pi;

Tx = cos(2.*pi.*Fsig.*t);               % Transmitted Signal
Carrier = cos(600.*t);                  % Carrier Signal
Modulated = Tx.*Carrier;                % Modulated Signal
Demodulated = Modulated.*Carrier;       % Demodulated Signal

Rx= filter([1 1],[1 -0.8],Demodulated);

figure(6);
plot(t,Tx,'r',t,Modulated,'b',t,Demodulated,'g',t,Rx,'k','LineWidth',1.5); grid on;
legend('Trans.','Modulated','Demodulated','Baseband')

How can I select the coefficients of the filter to perfectly recover my baseband signal?
Is there any other technique/command that can be used in matlab to recover the baseband signal other than this?


texts: What is the input dimension for LMS update algorithm

texts:
I just realized this forum existed, posted my question on the wrong forum.
Sorry about the formatting, please don't hesitate to ask if any part of this post is unclear. This is my first post.
I'm trying to implement a simulation of an ANC system with python, using this model here, but that's not the main point. My simulation keeps diverging, so I'm trying to figure out if it's my system is unstable or did I implement the NLMS algorithm wrong. Below is my code for the NLMS update. I'm using a source from Mathworks here.
x is input vector--------------------------------------------dimension(N * 1)
error is (desired vector - NLMS prediction)---------dimension(N * 1)
w_prev is the previous iteration's weight vector---dimension(1 * m)
def lms_update(x,error,w_prev,mu):

    ep = 0

    w_prev = np.array(w_prev)
    error = np.array(error)

    x_normal = np.dot(np.conjugate(x),np.linalg.inv((ep+np.dot(np.transpose(np.conjugate(x)),x))))

    w_update = w_prev + 2*mu*np.dot(error,x_normal)
    w_update = w_update.tolist()

    return w_update

then from above dimensions, after applying the update formula,
w(n)=αw(n−1)+f(u(n),e(n),μ)  and  f(u(n),e(n),μ)=μe(n)u'(n)
the update would be a 1 * 1 constant, and the weight would all be updated with the same value, which doesn't make too much sense to me.
So I adjust the input vector x into a (N * m) matrix, using vectors from current and previous iterations, so I can get the update to be a (1 * m) matrix. But I'm not sure if this is the correct way. From what I gathered from papers and google searches, nothing pops up indicating dimensions of the input. But I did implement a gradient descent algorithm before which input is a matrix, not a vector.


texts: Choosing Drive Signals for System Identification?

texts:
(I'm just learning a little about system identification so apologies in advance if this question is badly worded)

How do you go about choosing drive signals for system identification? I've seen PRBS signals used but it seems like that's going to work well for frequencies around the chip rate but not really low frequencies; I've also seen frequency sweeps.

If I have a SISO system that I know is close to a 2nd order linear system with poles in a certain range, and I can drive it with an arbitrary signal up to some amplitude A for up to some time length T, how do I pick a signal that would give me the best responses for determining the accuracy of the transfer function?

I tried googling for "system identification drive signals" but I don't see anything that pertains to my question.

edit: one particular type of SISO system I've dealt with is an (input=power dissipation, output=temperature) system for power semiconductor thermal behavior, and it seems very hard to model because there's usually a dominant pole at very low frequencies (<1Hz) and the next one might be 100 times higher, so any high-frequency drive signals just get very heavily attenuated.


texts: Why does the addition of noise improve coherence between two signals?

texts:
Taking a simple example of the coherence between two signals, each composed of two sinusoids of frequency 100 and 200 Hz, I expect the coherence to show two peaks in the spectrum located at precisely these two frequencies. The script below can be used to study the effect of adding noise An =0.1 or no noise An = 0. before taking the coherence. This produces in the case of added noise

and in the case of no noise

It is clear that with no noise there are many false peaks in coherence.

Question: Can someone show me mathematically speaking why this is the case?

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

t0, t1 = 0, 2.
fs = 2000
dt = 1./fs
t = np.arange(t0,t1-dt,dt)

#f1, f2 = 50, 50
#x1 = np.sin(2*np.pi*f1*t)
#x2 = np.sin(2*np.pi*f2*t+2*np.pi*0.)

An = 0.1            # Noise amplitude
f1, f2 = 100, 200   # Signal frequencies
p1, p2 = 0.25, 0.5  # Phase shifts (times pi)

x1 = np.cos(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + An*np.random.randn(len(t));
x2 = 0.5*np.cos(2*np.pi*f1*t-np.pi*p1) + 0.35*np.sin(2*np.pi*f2*t-np.pi*p2) + An*np.random.randn(len(t));

fcoh, coh = signal.coherence(x1, x2, fs, window='hamming', nperseg=400, noverlap=0)

fig2, ax2 = plt.subplots(1,1)
ax2.plot(fcoh, coh**2,"k-")
ax2.set_ylim([-0.1,1.1])
ax2.set_xlabel(r"$f$ [Hz]")
ax2.set_ylabel(r"Coherence")
ax2.grid("on")


texts:
What is LCCDE?I only know its abbreviation/full form :linear constant-coefficient difference equation

I know that in s domain we have differential equations and in z domain we have difference equation but what is LCCDE?What is its need?What is difference between a normal difference equation and a LCCDE?


texts: Aquila DSP C++ Library - Wave file FFT analog frequency off by factor of 4?

texts:
I'm working with the Aquila C++ DSP library. I'm computing the FFT of a wave file (16 bit depth, single channel, 44100 sample rate). I am using a window size of 16384 to calculate the FFT spectrum.

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

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

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

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

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

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

/* sample window */
int window_size = 16384;

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

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

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

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


texts: Implementing a Matlab filter to cancel cross talk between two sensors

texts:
I have two sensors that have a cross talk between them. I would like to cancel the cross talk. For this I recorded two tracks, where sensor no.1 (called x1) has some input, and sensor no.2 (called x2) is zero.

I was told to make the following operations:

  Since your cross talk is small:

  1) Simply measure the transfer function $H_{ba}(\omega)$ directly.

  2) Subtract it out as follows:

  2.1) Measure the transfer function from signal A to sensor B when signal B is 0.

  2.2) Create a filter from this transfer function (FIR or IIR, depending on it's shape).

  2.3) Now you can measure and subtract a filtered version from sensor signal A from sensor signal B:

  $y_b'(t) = y_b(t)-h_{ab}(t)*y_a(t)$. Where $h_{ab}(t)$ is the impulse
  response of your cross talk filter and $*$ the convolution operator.

Here is a Matlab code that I wrote for this purpose, but its performance is really bad. If someone has a suggestion what I did wrong, or can suggest a new code I will be grateful. If I may say, a friend helped me to implement this, and I think that we did not implement the filter part, but maybe I am wrong.

% The files designated to measure the transfer function.
% x1 has some input. x2 is is zero, and affected by the x1.

x1 = wavread(file1);
x2 = wavread(file2);

%The transfer function.

Hab = fft(x2)./fft(x1);

hab=ifft(Hab);

% Some new tracks, where x3 is recorded from sensor no.1 and x4 from sensor no.2.
% All the files have the same length.

x3 = wavread(file3);
x4 = wavread(file4);

c=conv(hab,x3);

% Subtracted signal

x4 = x4 - c(1:length(x4));

The following code did not subtract the cross talk on x4.
Here are some images explaining the problem. For some frequencies the above code actually did a good job and subtracted the cross talk, and for some it did the opposite and increased the energy. plot(abs(fft(x4))):

Good - Frequency energy subtracted:

Bad - Frequency energy increased (!):

Thank you!


texts:
To process a .wav audio file with numpy (using fast Fourier transform algorithm). I want to process an audio signal at a particular interval with a sampling frequency 44100hz and sampling rate of 20ms and also want to store the data into a array by using numpy with fft of 1024

Can anyone help me with the code?


texts:
I'm looking for a way to pan the PCM audio I have.
At the moment the audio is in a stereo format and I would like to pan it to the desired side (or more/less to a side).

I've recently figured out how to do mixing here so I guess it would be something similar (adding stuff to each other and such). What I was thinking is to add the audio data from one side to the other side.
For example if I were to pan hard left, I'd get the data from the right channel and add it to the left channel. I would then set the right channel to a constant value so that it would be silent.

Am I thinking in the right way or am I missing something obvious?

Thanks in advance!