Speech Recognition - Project Idea Wondering if someone could offer me any advice or help. Sorry if this question is too broad or anything.

I'm wanting to complete an algorithm in basic Speech analysis. The algorithm will detect whether an inputted sample is either saying yes, or, saying no.

I know that this has already been done, and I got my inspiration for this project from someone who has completed it using Zero-Crossing and in MatLab. His algorithm was as follows:

Divide the sound wave into evenly spaced blocks
Process each block for important characteristics, such as strength across various frequency ranges, number of zero crossings, and total energy.
Using this charateristic vector, attempt to associate each block with a phone, which is the most basic unit of speech, producing a string of phones.
Find the word whose model is the most likely match to the string of phones which was produced.

Now this algorithm / approach seems idealistic. But, the first point kind of confuses me. Ok, if my data sample is 4410 then how could I possibly split that into equal size blocks? For example, if the data size was 262144 then I know that multiples of 2 goes into it.

The second question is I want to use a FFT/linear predictive coding over zero-crossing but the FFT algorithm is really complicated to implement, so, would it be better to use a already built library?

My last question is this:

Do you think this project is possible, and, in under a year? I'm not completing it in MatLab, but, in C++. Also, will the above algorithm produce the results that I need? Any help would be really appreciated.

A few notes regarding your approach and detailed questions:

First, it is very common in audio analysis to split the signal to be analyzed into short overlapping windows, because audio signals are not stationary (their characteristics change over time), so they need to be processed on short segments over which they can be considered stationary (many analyses assume stationarity). Typical segment size is in the 5 .. 20 ms range. You'll be left with an odd-size segment at the end, just pad it with zeros and that's it.
In speech recognition the most widely features are MFCC. AR coefficients are quite common in older literature too. Zero crossing rate or energy are unlikely to discriminate phonemes in a robust way.

Now, regarding your project...

First, it is not clear if the goal of your project is to arrive at a working system even if you have to reuse existing code for that, or to do it yourself to learn on the road...

If you want to reach the goal without learning anything, there are ready to use speech recognition packages like Sphinx or HTK - the catch is that to understand these tools one actually has to know a lot about speech recognition tech. The setup time (training plans, model topology, grammar) and the concepts can be daunting.

An intermediate level would be to use third party code for signal analysis and feature extraction, and then focus on the recognition itself. You'll find C/C++ MFCC or AR modeling code in audio feature extraction packages like YAAFE or Aubio. You can also use the feature extraction command line utility of HTK (HCopy with the right config file for extraction + HList to get a clear text file).

When designing a speech recognition system, you have to ask yourself these questions:

What is the vocabulary size (from discriminating 2 words to recognizing natural language)?
What is the speaker variability (will I have to recognize from the same speaker again and again, or will I have to deal with many different speakers)?
What kind of training data will I have, in particular, is it possible to collect labelled speech samples from the same voice as the one you'll have to recognize?
What are the conditions (clear recordings... or voice through a glitchy mobile phone with background street noises)?

You have only answered the first question. If the answer to the other ones are these ideal conditions: the system will only deal with one speaker + will have a corpus of a few recordings of yes/no uttered by this same speaker + there will be no background noise ; the problem is not that hard. And these ideals conditions are not that unrealistic - there are some voice control PC applications which require this kind of prior training, the user needs to pronounce a few occurrences of each command word.

In these ideal conditions, a rustic approach (early 80s tech), but which will work, is DTW. Extract from each utterance in the training corpus (which will be labelled either "yes" or "no") a sequence of feature vectors, and compute the DTW distance between the utterance you need to recognize and those in the training corpus. Use this distance metric and a "k-nearest neighbour" scheme to decide on the label. DTW is archaic compared to a proper acoustic model + word model + language model graph, but for a very simple vocabulary (just a handful of words) it does the job.

Now, if you want something speaker independent (which will work for many speakers, including those for which we don't have a template), you'll have three plans. One is to stay with the same approach and record utterances from a very diverse set of speakers in terms of gender, accent, hoping that it will cover enough ground. But DTW scales very poorly! Another is to use an adaptation scheme - find the transform of your feature vectors that minimize the sum of DTW distance (or the minimum DTW distance) between your utterance and the database (paradoxically, I'm sure you'll find more literature about adaptation in the context of classic HMM-based recognition instead of DTW). The last approach is to forego DTW and do the real thing (phone-level GMM + word + language model).

I won't be as negative as the other members. I used to give a very similar project (digit recognition) to students which had no prior experience with audio (but a semester of signal processing and prior experience with Matlab + some background literature before starting the project), and they completed it in Matlab in two or three 3 hours lab sessions (roughly one for MFCC computation, one for DTW, one for the whole corpus management and k-NN + evaluation routines) with about 100 lines of code in total - so this is not an unimaginable engineering feat requiring a full team of PhDs...


Why does my signal magnitude increase after adaptive filtering? I am using a series-cascade of multiple NLMS adaptive filters each with step size 0.0040, leakage factor 1.0, and 100 filter taps. My signal gains magnitude at each step of the filtering, say the peak magnitude increases from 0.2 originally to 2.5 after using the first adaptive filter to 12.5 after the using the second adaptive filter on the output of the first one and finally 30 after using the third adaptive filter on the output of the second one. Why is this happening?

I have tweaked the step size and leakage factor but that did not help. Isn't my filter supposed to reduce the magnitude of my signal?

I don’t have the specific details for your filter but with digital filters in general it is typical for the filter to grow the signal in band in contrast to analog filters that shrink the signal out of band. It is all just a matter of scaling. Consider the simple case of a moving average FIR filter consisting of the summation of the previous N samples; such a filter will grow any low frequency signal within its passband, specifically growing a DC component by a factor of N. (We could then divide the result by N to get the actual average of try samples, but the primary filter itself here is growing the signal).

This approach of considering digital filters as "growing a signal in band" is particularly important for noise considerations in fixed point design, where you should typically avoid scaling the signal prior to the filter or the filters coefficients in order to normalize the result, but always allow the filter to grow the signal and then scale afterward. This is the reason for extended precision accumulators, and the reason is rather straightforward: if you scale the signal prior to the filter in a fixed point design you are effectively increasing the quantization noise. This quantization noise is typically modeled as independent white noise from sample to sample so in the filter will increase at a rate of $10Log_{10}(N)$ for the $N$ taps in the filter. If you scale after you are only adding this higher quantization noise level once.  Scaling the coefficients can be shown to accumulate quantization noise contributions in a similar way. This is also a reason why equalizers and adaptive filters should not be any longer than they need to be (to avoid noise enhancement).


Mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency? What is the mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency used?
This is partially related to the question that I asked on phase noise, but it seems that the complex magnitude of this component is the same for all PM signals using that carrier, just that the phase is rotated. Is it something like the average instantaneous phase difference?

What is the mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency used?

The phase of the carrier wave is your data-carrying signal. It also is the phase in complex equivalent baseband.
So, that's the mathematical relationship: equality.
Magnitude doesn't matter to PM. It might be the same for all signals you consider, but it doesn't have to be; it's simply irrelevant to the information signal.


External frequency divider in a PLL Does a frequency divider external to the loop BW of a PLL cause any phase ambiguity ? and why? and is the result of dividing by an odd number different than dividing by an even number?

Since you stated "external to the loop BW" I assume you mean specifically an output signal, not the signal fed back to the phase detector but provided to an external sink. Then yes, in this case a frequency divider would likely have phase ambiguity regardless of the division number unless special provisions are made to synchronize state (which requires an external time reference). This is because the division would start on any arbitrary edge (in typical cases the rising edge, but the salient point is any edge). Without an external synchronization of state, two identical PLL's may start the divided signal on different edges resulting in complete phase ambiguity. To eliminate this, some sort of state synchronization or reset is required.

Consider the operation of the divide by 2 as shown in the figure below. This typically operates by having the output change state on every rising edge of the input. Depending on the starting state of the local circuit, either of the two lower plots in the figure could result. Unless the circuit has a specific reference to common time there would be no way for it to remain in phase with its prior phase reference (were it to stop and start again, or considering two PLL circuits etc). This is what is meant by phase ambiguity. From this it should be clear that this would apply to any number N.


correct way for Zero padding I am doing simulation of FMCW radar. I am applying signal processing chain on raw radar data. Here in the code below, I am taking fft of data and plotting it. Moreover, I am doing zero padding by multiplying the signal length by 1000. I should be getting beat frequency at 2.020 MHz instead I am getting at 2.019MHz. There is something little probelem but I do not know what.

Is it the correct way for Zero padding? Moreover, am I defining frequency axes correctly? Thanks everyone in advance

% second-order butterworth
[b, a] = butter(2, [w1 w2], 'bandpass');

load('fb2030');                                   % loading the data
x = fb2030(30,:);
% filtering
y_filt = filter(b,a,x);                           % filtering the received signal
nfft = length(y_filt)*1000;                       %%%% Zero Padding
res = fft(y_filt,nfft)/ nfft;                     % normalizing the fft

f = fs/2*linspace(0,1,nfft/2+1);                  % choosing correct frequency axes
res = res(1:nfft/2+1);                            % amplitude of fft(taking the half length of nfft)
figure,plot(f,2*abs(res));
xlabel('Frequency in MHz');ylabel('amp')

This line seems suspicious to me:

f = fs/2*linspace(0,1,nfft/2+1);

Your frequency bins should be:

$$f_k = \frac{f_s}{\tt nfft} k$$

so

$$ \begin{array}{rcl}
f_{{\tt nfft}/2 + 1} &=& \frac{f_s}{\tt nfft} ({\tt nfft}/2 + 1)\\
&=& \frac{f_s}{2} + \frac{f_s}{\tt nfft}
\end{array}$$

whereas your linspace setting has it mapping to just $\frac{f_s}{2}$.


What are reds and trans in Image processing I'm trying to understand this code, however I'm unable to deduce the meaning of reds and trans. I would be glad if you can help with it:

import sys
import Image

# pixel constants in hex
trans="\x00\x00\x00\x00"
# reds most to least, the last is trans and used as a noise factor
reds=("\xf1\x34\x00\xff" , "\xfa\x00\x2f\xff",
   "\xf1\x36\x36\xff", "\xfe\x74\x74\xff", "\xfe\xc3\xc3\xff",
   "\xf9\xcb\xcb\xff", "\xf1\xd3\xd3\xff", "\x00\x00\x00\x00")

The complete code is at: https://matsu-project.googlecode.com/svn-history/r28/trunk/streaming-image-processing/src/reduce.py

The numbers appear to be RGBA values: red, green, blue, and alpha (transparency).

reds appear to be mostly red-tinted colors (except the last entry, which is a completely transparent  item).

trans is the transparent "color".


Finding delay parameters for watermarked signal I have a watermarked signal WATERMARKED SIGNAL

And the actual watermark WATERMARK

In Matlab, I need to find out the delay parameters $d_i$ which can be considered as the embedded
information.such that $0\leq d_i \leq 100$. There are 5 of them. I am afraid I don't understand the task, and my professor is on leave and will not be back before the task deadline. Would anyone be kind enough to at least tell me something about it? Show me a piece of code or a source? I can't seem to find any information about the code online.

Thanks.

All instructions:

EDIT:
Thank you Mathew for the answer. Now I believe I know what to do mathematically but still can't achieve the answer I am required (This must come from the fact I am an amateur with Matlab). I tried to play with your code but I have no idea what I need to change in order for it to work. So I will show you the code I created and if you are still interested in helping me, could you tell me what else needs to be changed?

clear
%Step 1. Upload sound files

[x,fs] = audioread('sample4_va18535.wav');
[m,fs] = audioread('assignment4_embedded_signal.wav');

%Step 2. Determine number of samples and create the window

N      = 132300;
w      = blackman(N);

%Step 3. Multiply both signals with the window to get x_w[n] and m_w[n]

xw     = x .* w;
mw     = m .* w;

%Step 4. Apply FFT to get X_w(k) and M_w(k)

XW     = fft(xw);
MW     = fft(mw);

%Step 5. Conjugate M_w(k) to get M_w*(k)

MWconj    = conj(MW);

%Step 6. Normalize XW * M_w*

Norm   = abs(XW .* MWconj);

%Step 7. Take IFFT

xwmw   = ifft((XW .* MWconj) ./ (Norm));

%Step 8. Plot everything (Using your method)

lags   = linspace(-N/2, N/2, N);

figure(1); hold off;
plot(lags,fftshift(abs(xwmw)))
% find the 5 peaks, then plot them
[pks,locs] = findpeaks(fftshift(abs(xwmw)),'NPeaks',5,'MinPeakHeight',0.0075);
hold on;
plot(lags(locs),pks,'rx');
xlabel('Correlation Lags'); ylabel('Scaled Correlation Magnitude (dB)')
grid on;
title('Watermark Cross-Correlation');

What I get is this : Surely it looks nothing like what you've got

Thank you for being patient with me. And I will understand if you can't find time to help me more.

Thank you for posting the additional information, it certainly helped.  I'll try to do a prose explanation with less "math" in it, and then point you in what is the right direction.

That's certainly a long winded problem statement, and I believe the solution is actually simpler that your professor may be making it out to be.  From reading the problem statement, a watermark has been made and was supplied to you.  This watermark is just sampled Gaussian noise, such that it lacks coherent frequency content, i.e. the spectral energy of the watermark is not located at any specific frequency (as opposed to say, a constant frequency signal like a sinusoid).

So we've got our original signal, we have our watermark signal, now it's time to actually watermark.  The simplest approach would just be adding the watermark signal to the original signal.  Your professor has gone a step ahead and added 5 copies of the signal.  Each one of these copies has its own time delay, as well as its own amplitude.  So, after we  add in our 5 watermarks to our signal, we have the signal Xm

Since there are 5 copies of the signal at 5 unique points in time in Xm, that means if we did a cross-correlation with our watermark signal and Xm, we would expect to see 5 distinct peaks at different amplitudes. This operation is typically called a matched filter, and you can think of it as using a finding the "locations" (in your professor's words, the delay parameters; think of this as when each watermark "starts" in Wm) of some given template (i.e. the watermark) in your recorded signal.

Your professor has gone ahead and detailed out how you can accomplish cross-correlation in the frequency domain.  This is often the preferred method, since it is typically much more computationally efficient.  Your professor has also detailed common ways to improve cross-correlation detection, specifically the use of windowing functions as well as the use of frequency normalization.

In turns out  those methods he described will be very important for your success.  The watermark signals are pretty low gain (naturally), so we'll need to maximize the signal to noise ratio (SNR) as much as possible.

I think that since this is seemingly an assignment of some sort, I'm not willing to just give you the answer in numeric or code form, but here's a pictorial output of what I got:

You can see that I've marked the 5 points corresponding to the time delayed version of the signal with red X marks.  These  peaks are easily 5x larger than any other point in the output of the cross-correlation, so one can easily  conclude that these are statistically significant, and point to the location of the delayed watermark signals.

Here's some code I used.  Feel free to play around with it and see if you get satisfactory results (note: "data" is your Xm signal, and "wm" is your watermark signal:

% Determine number of samples for the correlation lag vector
nsamps = numel(data);

% define a window function. Taylor windows give great performance
win    = taylorwin(nsamps,7,-10);

% window each signal
xin    = data;
wmw    = win .* wm;

% make the signals complex for phase coherency
xin    = hilbert(data);
wmw    = hilbert(wmw);

% perform cross-correlation in the frequency domain
xf     = ( fft((xin)) .* conj(fft((wmw))) );
% transform out of the frequency domain while normalizing
x      = ifft(xf ./ (abs(xf)+1));

lags   = linspace(-nsamps/2, nsamps/2, nsamps);

figure(1); hold off;
plot(lags,fftshift(abs(x)))
% find the 5 peaks, then plot them
[pks,locs] = findpeaks(fftshift(abs(x)),'NPeaks',5,'MinPeakHeight',0.0075);
hold on;
plot(lags(locs),pks,'rx');
xlabel('Correlation Lags'); ylabel('Scaled Correlation Magnitude (dB)')
grid on;
title('Watermark Cross-Correlation');

So hopefully that gets you going in the right direction.  Like I mentioned before, the watermark signals are pretty low gain, so make sure you do that normalization.  The windowing function is arguably not necessary, but they're pretty standard in practice since  they help your cross-correlation "behave" a bit better.


Expected value of a Delta function times a function I have a output signal $Y(t)= \sum_{n=-\infty}^{\infty} Z_n \delta (t-n\tau) h(t)$ and $Z_n$ is a random variable with equal probability to be $+-1$ and h(t) is $sinc(\frac t {\tau})$

and the Expected value of $Y(t)$ is

$E(Y(t)) = E[\sum_{n=-\infty}^{\infty} Z_n \delta (t-n\tau) h(t)]$

$=\sum_{n=-\infty}^{\infty}h(n\tau) E[Z_n]$

$=0$

But if I calculated in other way

$\int_{-\infty}^\infty t Y(t) dt$

$=\int_{-\infty}^\infty t \sum_{n=-\infty}^{\infty} Z_n \delta (t-n\tau) h(t) dt$

$= \int_{-\infty}^\infty  \sum_{n=-\infty}^{\infty} tZ_n \delta (t-n\tau) h(t) dt$

$= \int_{-\infty}^\infty  \sum_{n=-\infty}^{\infty} n\tau Z_n h(n\tau) dt$

$=  \sum_{n=-\infty}^{\infty} [n\tau h(n\tau)\int_{-\infty}^\infty   Z_n  dt]$

Why I have two different result with same signal? Thanks for any help.

If $h(t)$ is continuous at $a$, then $h(t)\delta(t-a) = h(a)\delta(t-a)$. Since the
$\operatorname{sinc}$ function is
continuous everywhere,
$$\begin{align}
Y(t) &= \sum_{n=-\infty}^{\infty} Z_n\delta(t-n\tau)h(t)\\
&= \sum_{n=-\infty}^{\infty} Z_n\delta(t-n\tau)h(n\tau)\\
&= \sum_{n=-\infty}^{\infty} Z_n\delta(t-n\tau)\operatorname{sinc}(n)\\
&= Z_0\delta(t)
\end{align}$$
since $\operatorname{sinc}(n) = \begin{cases}1,&n = 0,\\0,& n\neq 0.\end{cases}$
You (or your textbook's authors) probably meant to write something else.
You sure you are not missing a $\star$ indicating a convolution between
that $\delta(t-n\tau)$ and $h(t)$?


Does Short-time Fourier transform impact quality of signal? I am playing with spleeter to separate voice from an audio signal.
They go to the frequency domain by computing the Short-time Fourier transform to create a spectrogram. As there are some approximation (windowing and overlap of window) in this calculus, I was wondering if it hurts the quality of the signal when one computes stft and then inverse-stft.

STFT (short-time Fourier Transform) should not impact the quality of the reconstructed signal if properly implemented.
Furthermore, any signal $x[n]$ should be exactly recovered from its STFT analysis coefficients by a proper implementation of the associated inverse STFT;i.e., The forward and inverse pairs should exactly undo each other.
One should keep in mind about the overlap parameter and window functions however;i.e., if the window includes zeros, the inverse may not be possible at some overlaps. If the window does not contain any zeros, then the inverse is always computed.


the relationship between bandwidth and the datarate One Question of my note is like this.

Consider a communication channel with a bandwidth of 2400Hz. If QPSK technique is
used, what is the possible transmission rate? (Assume that rectangular
pulses are used in the baseband signals, and that 90% energy preservation is
required.)

For a baseband system, a bipolar signal of bandwidth B can support a rate of B
with 90% energy preservation (since bandwidth =1/τ in this case).
For a bandpass system, a BPSK signal of bandwidth of B can support a rate of B/2.
This is because, after modulation, the bandwidth is 2/τ.
For QPSK, rate can be doubled using the same bandwidth. Therefore the
supported rate is 2400bps.

Why the BPSK signal of bandwidth 2/τ can support a datarate of B/2? and what is the relationship between the bandwidth and the datarate? and why in the QPSK, the rate can be doubled using the same bandwidth?

Ideal symbol pulses are sequences of rectangular pulses.  Rectangular pulses, unfortunately, have horrible frequency profiles, as shown below.

The main lobe is where most of the pulse power is, and the other lobes are called side lobes.  Though they do diminish in power the farther out they go, they don't diminish very quickly.  The main purpose of pulse shaping is to dramatically reduce the power of those side lobes.

The width of the primary lobe is $\frac{2}{T}$, where $T$ is the rectangular pulse width, which is equivalent to the symbol period.  Thus, once pulse shaping is employed, the signal bandwidth is $\frac{2}{T}$.

Since BPSK transmits one bit per symbol, and the symbol rate is $\frac{1}{T}$, BPSK can transmit $\frac{1}{T}$ bits per second.  Therefore, the bit rate divided by the bandwidth is $\frac{\frac{1}{T}}{\frac{2}{T}}$, which reduces to $\frac{1}{2}$.  Thus, the bit rate is $\frac{1}{2}$ the bandwidth.

When the signal is at baseband half of the mainlobe is in the positive frequencies and half is in the negative frequencies.  Given that it is a real signal, the two halves mirror each other so some people consider it the bandwidth to be halved, or $\frac{1}{T}$.  Thus, at baseband the bit rate would be equal to the bandwidth.

QPSK transmits two bits per symbol, so the bit rate for QPSK is $\frac{2}{T}$.  It follows that QPSK can transmit $2$ bits per Hz of bandwidth at baseband, and $1$ bit per Hz at passband.


how sobel edge detector as a first order derivative is turned to a 2 dim filter? I want to know theorem behind sobel operator. basically a first order derivative is -1/2 0 1/2 . its a vector. the question is how this vector is turned to a matrix below and what was theory and reason behind that?
-1 0 1;
-2 0 2;
-1 0 1;

Background

The Sobel edge detector was introduced back in 1968 by Irwin Sobel and Gary Feldman as the Sobel-Feldman operator. In broad strokes, 'edges' in images are related to gradients, which motivated their development of a discrete differentiation operator. The Sobel-Feldman operator only computes an approximation of the gradient and not the actual gradient. However, this integer valued operator is simple and fast, which makes it useful.

  Directionality

The Sobel-Feldman operator is only capable of estimating gradients along mutually perpendicular directions separately and not jointly. Since an image is typically a matrix itself, the operator provides gradients along rows and gradients along columns. This limitation is important since we will use this property to extend the operator to multiple dimensions.

  Matrix Decomposition

The Sobel-Feldman operator in 2 dimensions is given as a kernel matrix. For image matrix $\mathcal{I}$, the Gradient operator (Sobel edge detector) for x- and y- directions is given as:

$\mathbf {G} _{x}={\begin{bmatrix}+1&0&-1\\+2&0&-2\\+1&0&-1\end{bmatrix}}*\mathbf {\mathcal{I}} \quad {\mbox{and}}\quad \mathbf {G} _{y}={\begin{bmatrix}+1&+2&+1\\0&0&0\\-1&-2&-1\end{bmatrix}}*\mathbf {\mathcal{I}} $

Sobel kernels can be decomposed as the products of an averaging and a differentiation kernel. For example,

$\mathbf {G} _{x}= {\begin{bmatrix}+1&0&-1\\+2&0&-2\\+1&0&-1\end{bmatrix}}={\begin{bmatrix}1\\2\\1\end{bmatrix}}{\begin{bmatrix}+1&0&-1\end{bmatrix}}$

This is the key idea in extending the Sobel operator from 1-D vector to 2-D matrix. The operator uses a triangle-kernel for averaging (smoothing) in the direction perpendicular to the desired gradient direction and the Sobel differentiator in the desired direction.

  Higher Dimensions

Sobel–Feldman filters for image derivatives in different dimensions with $x,y,z,t\in \left\{0,-1,1\right\}$:

vector: $h_{x}'(x)=h'(x); h_{x}'(x)=h'(x)$

matrix: $h_{x}'(x,y)=h'(x)h(y) h_{x}'(x,y)=h'(x)h(y)$

Tensor: $h_{x}'(x,y,z)=h'(x)h(y)h(z) h_{x}'(x,y,z)=h'(x)h(y)h(z)$

Tensor:  $h_{x}'(x,y,z,t)=h'(x)h(y)h(z)h(t) h_{x}'(x,y,z,t)=h'(x)h(y)h(z)h(t)$

  Tensor Sobel operator

For 3-D (x,y,z), the Sobel-Feldman operator in z-direction is:

$h_{z}'(:,:,-1)={\begin{bmatrix}+1&+2&+1\\+2&+4&+2\\+1&+2&+1\end{bmatrix}}\quad h_{z}'(:,:,0)={\begin{bmatrix}0&0&0\\0&0&0\\0&0&0\end{bmatrix}}\quad h_{z}'(:,:,1)={\begin{bmatrix}-1&-2&-1\\-2&-4&-2\\-1&-2&-1\end{bmatrix}}$

  Summary

The Sobel operator is only a poor approximation of gradient using a center-difference operator.
The operator is only informative about 1 direction at a time.
The operator is extended to higher dimensions by using a triangle shaped averaging filter for other perpendicular directions (dimensions).
There isn't a deep theoretical reasoning behind this operator. It is a discretized version of a differential (gradient) and Gaussian filter (smoothing). But it is simple and fast.


Estimating harmonic loudness I'm working with data that looks something like this:

sec Hz amp
0.0 2630 0.001
0.1 2340 0.1
0.2 2220 1.0
0.3 2200 0.7
0.4 2200 0.2
0.5 2200 0.01

sec Hz amp
1.0 463 0.001
1.1 534 0.1
1.2 440 0.4
1.3 440 0.7
1.4 440 0.2
1.5 440 0.01

This data is for two separate harmonics. You can see that the frequency drifts a little with the attack of the note and that the amplitude beings to fade out towards the end of each harmonics. I'm not showing all the samples.

Here's what I'm doing to establish the relative loudness of these harmonic.

Harmonic 1 has a maximum amplitude in the frequency domain of 1.0.
Harmonic 2 has a maximum amplitude in the frequency domain of 0.7.
So there amplitude ratio is equivalent to 3.1 dB.

I now take into account that loudness is dependent on frequency and do so via the following table:

(freq,dbs)
(0, -50)
(20, -43)
(30, -33)
(40, -27)
(50, -23)
(60, -21)
(70, -19)
(80, -17)
(90, -16)
(100, -15)
(200, -8)
(300, -6)
(400, -6)
(600, -6)
(700, -7)
(800, -8)
(1500, -9)
(2000, -7)
(2200, -6)
(2500, -4)
(3000, -1.5)
(3700, 0)
(4000, -0.5)
(5000, -4)
(6000, -9)
(7000, -14)
(9000, -16)
(10000, -15)
(12000, -25)
(15000, -40)
(20000, -55)
(30000, -70)

Harmonic 1 has a frequency of 2200 Hz, which gives a value of -5 dB in the look up table.
Harmonic 2 has a frequency of 440 Hz, which gives a value of -6 dB in the look up table.
So that's a difference of 1 dB.

The total difference is then 3.1 dB + 1 dB = 4.1 dB. Right?

So a difference of 10 dB is equivalent to a doubling of loudness:

2^(dB/10)=increase in loudness, with
2^(4.1/10) = 1.33

So the first harmonic is approximately 33% louder than the second.

Am I on the right track with regard to estimating the difference in loudness? Any obvious mistakes. Any obvious improvements?

It would help if you could state what exactly you are trying to do.
The "relative loudness of the harmonics" is normally a function of time because of the different decay rates and/or temporal envelope of the harmonics. There are different ways to condense this into a single number: initial amplitude, peak amplitude, total energy, 250ms into the decay, etc. The most suitable method would depend on your specific goal.

The look up table looks questionable to me. None of the loudness curves that I'm aware of (A-weighting, C-weighting, Fletcher Munson, Robinson Dadson, ISO 226) look like that. -9dB at 1500 Hz seems just plain wrong.

The shape of the loudness curves is a function of absolute level. So whatever curve you pick it can only be correct at one level (or a small range of levels). In order to get this really right, you would need measurements that are calibrated in Pascal or dBSPL.


Phase of FMCW linear chirp I'm doing a project involving FMCW linear chirp signals.

Given the phase

$\phi = 2\pi\frac{k}{2}t^{2} + 2\pi f_{c} t + \theta$

where $k$ is the frequency increase rate.

When transmitting frequency modulated pulses continiously will the phase increase in this manner?

Or simply repeat in this manner?

What is the advantages and disadvantages of either one?

Remember that phase is usually only defined on the $[0;2\pi[$ interval (or on $[-\pi;+\pi[$). When you do that $\mod 2\pi$ operation on your graphs, you'll notice that there's no significant difference between your two plots.


How to generate a sinusoid with STFT at a frequency which is between bins Is it possible to generate a sinusoid at a frequency which is between two bins (so a fractional index) with the STFT algorithm (by filling bins and doing the inverse) ? If yes then how can this be done ?

In my experiment with overlap add STFT algorithm and with a window size of 1024 i just compute the distribution of energy between the two bins that the sinusoid fall into however this produce two sinusoids plus some noise in between... is that normal ?

A sinusoid with a frequency that is between bins in the FFT frequency domain is circularly discontinuous in the time domain.  So you can't use the same IFFT results back-to-back without the noise from this discontinuity between each IFFT window, as the end of one window will have a value too far from the beginning of the next window (unless your frequency is bin centered, e.g. exactly integer periodic in the FFT's length).

If you calculate and vary the starting phase for each successive (complex input) IFFT, then you might be able to reduce or get rid of noise from the window discontinuity.

Note that the spectrum of a "between bin" frequency sinusoid is not just a cross-fade between the two nearest bins, but a Sinc shaped spectrum, spanning the full width of the FFT.  (Or a windowed Sinc function for a decent approximation.)


Drawing the modulus from a Transfer function As you may guess from my other questions i am a student so pardon me any ignorance and guide me to the truth. Thanks

In this exercise from an exam i was given that transfer function with poles and zeros. I was asked to draw the modulus and the "UHD imagive below the Transfer function" image was my answer.

I know that when you encounter a zero, you "go down" and when you encounter a pole you "go up" in the graph. When you encounter both "aligned", you get a cero as they cancel each other.

The answer from the teacher was that the maximum, which is 1, is at pi/8 (and -pi/8), not at pi/4, because at pi/4 (and -pi/4) he draw the minimum.

Why is that? and what i did got wrong or missed in my understanding?

Thanks for answer.

EDIT: Here i add another example from a recent exam. I would like if someone could do a brief/short explanation in how to properly draw it and what should i know or notice.

The question is the same, but i lack proper understanding of it.

Your sketch is not bad at all. I don't know why your teacher drew a "minimum" (do you mean a zero?) at $\pi/4$, because there shouldn't be any. The zeros are clearly at $\pm\pi/2$ and at $\pi$, and the main peak should indeed be close to $\pi/4$. By the way, I'm missing a pole. There are $4$ zeros and only $3$ poles. The fourth one should probably be at the origin. Maybe your teacher didn't give this exercise as much thought as he should have ...

EDIT:

Concerning your second example, it's always the same: a pole gives a bump in the magnitude response, the closer the pole to the unit circle, the more pronounced (sharp) is the bump. A zero on the unit circle gives an exact zero in the magnitude response, a zero away from the circle gives a dip. The closer the zero to the circle, the more pronounced (deeper) is the dip. Note that the exact maxima and minima in the magnitude response caused by poles and zeros close to the circle are often not exactly at the frequencies that correspond to the angles of the poles or zeros, but often they are close. This is all about a basic sketch, so no need to try to be overly precise.


Why can you use the one-sided laplace transform to solve differential equation describing a causal LTI-system? In an example, an equation describing a causal LTI-system is
$$
(D^2 + 5D + 6) y(t) = (D+1) x(t)
$$
where $y(t) = y_{zs}(t) + y_{zi}(t)$ and the initial conditions are $y(0^-) = 2, \dot{y}(0^-) = 1$.
$x(t) = e^{-4t}u(t)$ and we want to calculate $y(t)$.
My teacher said in an example of a solution that because $x(t) = 0, \quad t < 0$ we can take the one-sided/unilateral laplace transform of the RHS (since then the unilateral and bilateral laplace transform are the same, I guess ). Further, the explanation continued with that because the system is causal, the impulse response $h(t) = 0$ for $t < 0$ and therefore $y_{zs}(t) = (x*h)(t) = 0$ for $t < 0$. Therefore, supposedly, $y(t<0)=0$ and we can take the one-sided laplace transform of the LHS too.
I am questioning the part in boldface, because what about the zero-input response $y_{zi}(t)$, how do we know its value for $t<0$? I would like to belive it is not equal to $0$ for $t<0$ due to the initial conditions given at $0^-$ not being zero and $0^- < 0$. If it is not equal to $0$ for $t<0$ how can we know we can take the one-sided laplace transform of the LHS?
Note: I also read this question in which it is stated that non-zero initial conditions make the system non-linear and time-varying, which also makes me think the equation at hand is inconsistent with the fact that it describes an LTI-system?
Edit 1 in response to the comment about Lathi's linear systems and signals:
I did not remember that this example was indeed from the book, but I have read the example and the section "Comments on initial conditions at $0^-$ and at $0^+$ " before. The section explains that we cannot expect the total response $y(t)$ to satisfy the initial conditions given at $0^-$ at $0$, which makes perfect sense I think since $0^- \neq 0$. It goes on to say that there is another version of the laplace transform, $L_+$ that is not as convenient to work with. The author's discussion might have gone a bit over my head, because unfortunately I can't understand how it answers my questions, that is how we can assume $y(t)$ to be causal (in order to use the unilateral laplace transform on both sides of the equation) and the note about the non-zero initial conditions implying the system to not be LTI.
The reason I have not accepted the answer given so far is that it is a bit to advanced for me to judge its correctness and therefore I wanted to wait a bit in case there would be more input for my question. But eventually I will just assume it is correct and accept it.

If it is not equal to 0 for t<0 how can we know we can take the one-sided laplace transform of the LHS?

Unfortunately, the way the problem is framed it is very difficult to do in a rigorous manner.  You can, by remembering that $0^-$ is shorthand for $-\epsilon$ as $\epsilon \to 0$, and taking limits everywhere, and having expressions where you're looking at intervals of $t$ from $-\infty$ to $-2\epsilon$ and other maddening things.
The easy way to do it is to separate the problem: solve the problem for the $x(t)$ that's given.  Then augment the problem by finding an $x(t)$ that is zero everywhere except in $0^- < t < 0^+$ that will result in $y(0^+) = 2$ (note the $0^+$ instead of $0^-$) and $\dot y(0^+) = 1$.  This involves setting $x(t)$ to a linear combination of $\delta(t)$ and $\delta^2(t)$ (as if $\delta(t)$ weren't wacky enough, $\delta^2(t) = d/dt\ \delta(t)$).
Then add the two solutions together.  It may not be what your prof had in mind, but it takes the least amount of hand-waving, and if you really want rigor you can break $\epsilon$ and put it to work.
By doing the above trick of expressing $x(t)$ as a linear combination of the specified $x(t)$ plus whatever sum of the Dirac impulse and it's derivatives, you cast the problem into one where the system is linear -- you're just using physically impossible values for $x(t)$.


Help in problem formulation for estimation of image as a feature vector - SISO or MIMO FIR channel model? Based on the paper Blind Image deconvolution:

A feature vector is a list of numbers used to represent an image. The feature vector for my case takes values as symbols $-1,1$.

An instance or an example of an image is represented as a feature vector of $d$ feature values $\mathbf{x} = [x_1,x_2,...,x_d]$. I have $N$ examples
$$\{\mathbf{x}_i\}_{i=1}^N$$ in the form of a database. So, the database contains $N$ rows of feature vectors each of $d$ length.  One such pattern (feature vector) is
$$\mathbf{x} =  [1, -1, -1, -1, -1, 1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, 1, -1, -1, -1, -1, 1]'$$
where $d = 30$
The database thus has $N$ rows of examples each of $d$ length :
$$ \mathbf{x_1} = x_{11},x_{12},....,x_{1d}$$
$$ \mathbf{x_2} = x_{21},x_{22},....,x_{2d}$$
$$ ::$$
$$ \mathbf{x_N} = x_{N1},x_{N2},....,x_{Nd}$$

I don't know how to use the above as input to the channel. In general,
the output is,
$$\mathbf{y} = H^T x$$

Assuming that the impulse response is that of FIR, then I am facing difficulty in how to represent the channel. Should the source input be each component of the feature vector or an entire example?

if the input in the form of the feature vector is transmitted via a channel whose impulse response is modeled as moving average (finite response, FIR) then inputs and how many outputs should there be? The output would be the convolution of the input and the impulse response. Would the channel be Single input-Single Output FIR or Multiple Input Multiple Output FIR ? The objective is to estimate the channel parameters and the input using estimation methods from the output observation only. Would the estimation be performed using one example or all the examples?
Should there be $d$ channel coefficients?
A mathematical representation of the system will really help to clear the concept. Thank you

This question is confusing, but I thought I'd try to make some sense of it.

Suppose you are transmitting one of your $\mathbf{x}_i$ vectors one component at a time. That is, the transmitted signal, $s[k]$,  for the $i^{\rm th}$ feature vector is:
$$
s_i[k] = x_{ik}
$$
If this signal is corrupted by a channel with impulse response $h[k]$ then the received signal, $y[k]$ will be
$$
y[k] = h[k] \star s_i[k] = \sum_{m=0}^{M-1} h[m] s_i[k-m] = \sum_{m=0}^{M-1} h[m] x_{im}
$$
where $\star$ represents convolution.

If we use this formulation, then the answer to:

  Should the source input be each component of the feature vector or an entire example?

is that it should be each component of the feature vector.

  if the input in the form of the feature vector is transmitted via a channel whose impulse response is modeled as moving average (finite response, FIR) then inputs and how many outputs should there be?

There should be a sequence of outputs. If the FIR filter is length $M$, then the output will be of length $M+d-1$.

  Would the channel be Single input-Single Output FIR or Multiple Input Multiple Output FIR ?

In the model above, it is a SISO system.

  Would the estimation be performed using one example or all the examples?

If the channel does not change, then you might be better off using all the data you have available.

  Should there be $d$ channel coefficients?

There is no point in tying the number of feature vector components with the length of the channel. In the above, the length of the channel is $M$.


Generating a signature to identify RF modulation types Thanks to the Realtek 2832U I'm just discovering the world of Software Defined Radio.  There are a lot of signals I'm able to pick up with this device, but identifying them is not always easy and cannot be done automatically.

When scanning for radio transmissions, I can look at a waterfall plot and see whether a broadcast is AM, FM or a handful of digital modes I have learnt to identify, but this is slow and error prone.  I am wondering whether there is a way to automate this.

What I am after is an algorithm to examine the block of RF spectrum that has been digitised, to identify any signals along with their transmission modes.  For example, in a particular part of the spectrum it might identify three broadcasts, two that are FM and one that is a POCSAG data stream.

I am thinking that for this you would need to generate some sort of signature for each detected signal, and then compare that to a list of known modulation types.  Is it possible to generate a signature like this, taking into account the different amplitudes (signal strength) that result from varying distances to the transmitter?  The signature (i.e. an identifying number) would need to be the same for all instances of that signal - for instance, each consumer FM radio station would have the same signature.

I am not really sure where to start, or how to deal with the bandwidth ranges involved (some signals are narrow and close together, others are very wide.)  I'm also unsure whether you would have to pick out strong signals first and only process those (losing weaker ones) or whether you would examine the spectrum uniformly regardless of the power level at that particular frequency, discarding any frequencies that matched the signature for "no signal present."

The reason I am asking is that I would like to write a program that can tune the radio to a given frequency and then display all the transmissions it picks up, decoding them all automatically where possible.  In order to do this, it must be able to identify signals in enough detail to demodulate them.

Yes, this is quite a studied problem in certain fields. I co-own a patent on it and we used one out of many approaches. We broke the problem into a feature extraction step and then a classification step. Some features can be the bandwidth, the continuity of the phase, the amount of phase clusters for digital modulation, etc. Then we used a Mahalanobis distance comparison to known modulations. I've also seen linear discriminants and neural networks used to process the features.

Certain modulations also exhibit certain shapes once you square them, cube them, etc. and then compute the power spectral density. Often if you do an FM or AM demodulation you can see spectral lines for digital modulations. Cyclostationary computations can also be helpful in digital modulation detection. The post above is right in regards to the search topic: "modulation classification"


What is the type of each 4 different noise signals that are shown below? Dear signal processing experts,

I am looking at four different noise signals.I got them in a time domain and they are different. I went to the frequency domain and found the fft for all four signals.I got differnet shapes but I still cannot recognize which type of the noise(Pink,blue,white brown or what type)[![enter image description here][5]][5]

Could you please help me? what should I do?

Let me start with Noise 4. Noise 4 is "Zero-mean White Noise".  White noise means the signal power is distributed across all frequencies evenly, which is what you see in the FFT result. From the histogram of the time signal we see that the signal is almost uniformly distributed, so close to a uniform white noise process. Although it appears closer to the sum of two independent Gaussian distributions, one with a positive non-zero mean and the other with a negative non-zero-mean, approximately +/- 0.25. (The result would be a zero-mean white noise random process).

Noise 1 is non-zero mean white noise, as evidenced by the spike at $F=0$ with a similar white noise spectrum everywhere else.  From the histogram we see that it appears to be the sum of a zero-mean Gaussian White Noise (AWGN) process with a non-zero mean random process. The power level in the non-zero mean process is so small compared to the zero-mean AWGN process that we are unable to really discern further properties of that non-zero part that is added except for the mean which appears as a strong component in the FFT. Also from the histogram it looks like there really is very little signal component between the two, so the time domain plot appears as it does simply because you are connecting the dots as a line plot rather than plotting just the points (try using plot(time, noise, '.') to see just the waveform points. I suspect that it is a Gaussian distribution with a mean of 10 that is sparse compared to the non-mean noise waveform.

Noise 2 has been low pass filtered, and looks like it may have possibly been a white noise process that has been filtered by adding to itself after having been delayed by one sample. This could be called Pink Noise since the low frequency noise dominates (drawing an analogy from the visible light spectrum). From the histogram we see that this is well approximated as a Gaussian random process. So in this case it would be band-limited Gaussian noise.

Finally Noise 3 has a negative DC offset. Because of the scale, we cannot see if the noise process is pink or white (would need to zoom in on the FFT to observe that.) However from the histogram we see that it is a Gaussian distribution with non-zero mean.


ICA for Noise Reduction over a Dataset Suppose my dataset consists of $N$ example vectors $\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}$ where $\mathbf{x}_{n} \in \mathbb{R}^{p}$ $\forall n$. I assume that each vector  $\mathbf{x}_{n}$ is comprised of an underlying true datapoint $\mathbf{s}_{n}$ that is corrupted by additive Gaussian noise $\mathbf{x}_{n} = \mathbf{s}_{n} + \mathbf{w}_{n}$.

What I want to do is estimate the $\mathbf{s}_{n}$. Can I do this with ICA?

So the ICA model is of the form:

$$ \mathbf{x}_{n} = a_{1,n} \mathbf{s}_{1} + a_{2,n} \mathbf{s}_{2} + \ldots + a_{N,n} \mathbf{s}_{N} $$

Does this mean I need to fit my problem to the model by assuming that $ \mathbf{s}_{n} = a_{1,n} \mathbf{s}_{1} $ and $\mathbf{w}_{n} = a_{2,n} \mathbf{s}_{2}$ ? If so, how do I recover the vectors $\mathbf{s}_{n}$ $\forall n$?

I ask because I see the ICA problem formulated as:

$$ \mathbf{x} = \mathbf{A} \mathbf{s}  $$

Would this then imply that: $\mathbf{s} = \mathbf{A}^{+} \mathbf{x}$, where $\mathbf{A}^{+}$ is the psuedo-inverse of $\mathbf{A}$? But how would I get the $\mathbf{s}_{n}$ out of this formulation?

Note: I cross-posted this on stats.stackexchange to get the statisticians' point of view, but I'd like to hear what the signal people have to say.

Edits:

So I do not know $\mathbf{A}$ but I do know that each $\mathbf{s}_{n} \sim f_{s}(\cdot | \theta)$ for some parameter vector $\theta$. I have specified $f_{s}(\cdot | \cdot)$ but $\theta$ is unknown. I also know that the noise vectors are not sampled from that same distribution, but that they are independent of the signal vectors.

This appears similar to a classic least squared solution of an overdetermined equation that proceeds as follows:

Starting with:

$$ \mathbf{x} = \mathbf{A} \mathbf{s}  $$

$\mathbf{A}$ is not a square matrix if overdetermined (more equations than unknownns) so therefore an inverse does not exist. What you do then is multiply both sides by the transpose of $\mathbf{A}$ since $\mathbf{A^TA}$ is a square matrix. From that, assuming and inverse exists for$\mathbf{A^TA}$, you can solve for s:

$$ \mathbf{A^Tx} = \mathbf{A^TA} \mathbf{s}  $$

$$ \mathbf{s} =  (\mathbf{A^TA})^{-1}\mathbf{A^Tx} $$


Computing the derivative from a set of real-time data? So I have an array of say 25 samples and I would want to be able to note the trends of whether it's decreasing n or increasing from those 25 sample time interval(basically 25 samples array is my buffer that is being filled by every say 1 ms). Note that it is general trend that I am looking for, not the individual derivative(as I would have obtained using finite difference or other numerical differentiation techniques). Basically I expect my data to be noisy so there might be ups and downs even after doing filtering and so on. But it's the general trend of increasing or decreasing behaviour that I am looking for. I want to integrate the increasing/decreasing behaviour in every ms to trigger some event which is more of a user interface event (blinking a LED) so it does not have to very delay of processing as long as I can detect the general trend.

Thanks in advance!

Depending on the level of noise you expect to encounter on your signal, you might actually be able to use the finite differencing you suggest, originally. If you want to know the general trend of your data upward or downward, sum together the, in your case, 24, differences. So, lets say that you have some metric, $D$, for your signal. Lets say that for a given signal, $x[n]$, the first order difference is denoted as $x\prime[n]$ (and is, of course, one entry shorter than $x[n]$). Define $D$ as

$$ D = \sum_{n=0}^{N-2} x\prime[n] .$$

Then, classify your signal as uptrending or downtrending via

$$ D > 0 ~\therefore~\uparrow$$
$$ D < 0 ~\therefore~\downarrow$$

I made some basic assumptions about the nature of the signals in question. I just did the most basic analysis and thought, well, maybe just use a ramp, $r[n]$, corrupted by some noise, $z[n] \sim \mathcal{N}(0,\sigma^2)$. This may or may not be the situation you're looking at, so you'll need to adjust accordingly and see what works for you. But my final model was just the simple
$$x[n] = r[n] + z[n]$$
for the length you gave, $N = 25$.

I went ahead and ran some experiments in MATLAB. I used the following code:

experiments = 1000;
N = 25;

slope_min = -2;
slope_max = 2;
slope_range = slope_max - slope_min;
noise_level = 2;

for i=1:experiments

    slope(i) = slope_range*rand(1)+slope_min;
    r = slope(i)*(1:N);
    z = noise_level.*randn(1,N);

    x = r + z;

    diffsum(i) = sum(diff(x));
end

First, let me show what some of the signals look like...

So, I varied the variance of the noise signal, $z[n]$ and recorded the following correlations between the true slope (up or down trend) of $r[n]$ in comparison with the metric we defined above, $D$. If $D$ does a good job determining the trend of the signal, we should see a strong correlation between downtrends (negative slopes) and negative $D$ and uptrends (positive slopes) and positive $D$.

We can see from the above that $D$ does seem to do a decent job. As the noise level increases (or the signal moves further away from our model), the correlations are not as strong. Try it out on your signals, see if does anything for you.


Performing FFT in Matlab set of signals containing time and amplitude was given to me and asked to process the signal with the help of FFT function. I need help in converting the signal to frequency domain and to extract the feature from the signal so that I can compare the test signal with the given signal so the signal is a fault signal or not

If doing this in Matlab the command is fft(x) where x is your signal. Note that the result will be presented as a complex frequency vector of N points (or bins) from 0 to N-1.  Bin 0 is representative of the DC term and bin N-1 is 1 bin less than your sampling rate. Due to aliasing in digital sampling, this also represents the first negative frequency bin. Often we wish to view the spectrum with DC in the center with the positive and negative half spectrums. To do this easily in Matlab use fftshift:

out = fftshift(fft(y))

Also note that out will be complex, so you may want to view magnitude and phase separately:

mag = abs(out)

or

magdB = 20*log10(abs(out))

phase = angle(out)    <- in radians.

Hope that helped!


How to Remove the Patch Artifacts of Neural Network Denoising Process? I have written a python script which uses the Noise2Noise: Learning Image Restoration without Clean Data implementation of the Auto Encoder which is useful to remove noise from images. In the original paper implementation they were using full images so there were no border artifacts.

I want to use my script on raw images from long exposures and high ISO to remove the noise. For this purpose, I fragment the image into small fractions which are run through the Auto Encoder network. After prediction, original image is reconstructed. However, the problem is with the border connections between slices (As showed in the images below). Do you have any suggestions what would be the best idea to "fix" it? I would be more than grateful!

One simple way to solve it is using Overlapping Patches.

Let's say you have image which is $ 20 \times 20 $ and you work on patches of the size $ 5 \times 5 $.
As I understand from your description you do 16 times denoising of $ 5 \times 5 $ patches.

What you should do is run the patches mask like in convolution. So each pixels (Ignoring boundaries) will be filtered in 25 different patches.

Now you have 25 different estimates for each pixel.
The trivial solution is to use their mean.
More advanced methods will use some prior on how to create a good weighted combination as output.
One prior could be keeping the composition of the image as smooth as possible.
Other methods could be built around estimation of the noise in each patch.


Why do we need to flip the kernel in 2D convolution? Why do we need to flip the kernel in 2D convolution in the first place?
What's the benefit of this?
So, why can't we leave it unflipped? What kind of terrible thing can happen if you don't flip it?

SEE:
"First, flip the kernel, which is the shaded box, in both horizontal and vertical direction"

http://www.songho.ca/dsp/convolution/convolution2d_example.html

At the risk of repeating what has been said elsewhere, this is how convolution in two dimensions is defined:

$$y(n_1,n_2) = \sum_{k}\sum_{l}x(k,l)h(n_1-k,n_2-l))$$

$x(n_1,n_2)$ is your data and $h(n_1,n_2)$ is the kernel. Note that in order to compute the double sum above, you need to flip $h(n_1,n_2)$ with respect to $x(n_1,n_2)$. Just take $n_1=0$ and $n_2=0$: in this case you need to multiply $x(k,l)$ with the flipped kernel $h(-k,-l)$. In order to compute the output for other values of $n_1$ and $n_2$ you need to move the flipped kernel across the input data.


Generating inaudible sound waves using MATLAB I want to generate inaudible sound waves (above 18kHz) with 48kHz sampling rate using MATLAB.

Generating a sound wave with 18kHz frequency makes it inaudible.

I've inserted the data in the 18kHz band in the frequency domain, converted it into time domain signal.
Then I've duplicated the signal, concatenated it, and played it.
As expected, the played sound was inaudible.

The problem is,
I want to play the inaudible sound with an interval.
(playing sound wave for 1.75ms and having a silence period for 4.17ms
 repetitively)

Because of the inserted silence period, the concatenated signal is no longer inaudible because new frequency components are added.

How I can transmit inaudible sound with time interval?

Eventhough the exact cause of those wide-spread line spectra is not very clear to me from the supplied information, it's most probably due to the on-off switching implied by the silence period you have added in between your message signals.

The on-off waveform is an implicit operation you multiply to your signal. Which has a value of $1$ during the sound interval and a value of $0$ during the silence interval. The fundamental period of that waveform is roughly 200 Hertz.

More technically, the on-off switching with a given duty-cycle will shift the message spectrum by the amount of k-th harmonic frequency for each significant harmonic present in the continuous Fourier series of the on-off waveform.

In fact this is a method of realizing an analog modulation algorithm. Multiplication with a rectengular wave creates a multitude of harmonic centered up-mixed spectrums, of which only one is selected by a suitable bandpass filter.

Therefore as a solution offer, you could try the following:

1- Generate the 18kHz (already bandpass) message signal in time domain.

2- Insert the silence periods accordingly.

3- Apply a BPF (bandpass filter) centered around 18Khz.

When implemented carefully this should yield a solution to your problem I guess. Distortions in the message signal can happen however.


How to estimate covariance matrix using Fourier representation? So, I have multidimensional time-series $X \in R^{(d \times T)}$, and I want to determine the covariance matrix of that signal in a specific frequency band.

I might filter the signal to that specific frequency band and do the usual $XX'$. However, if I assume the signal to be stationary, I might use the Fourier representation of the signal for a much efficient implementation, according to the Plancherel theorem

However, using only a subset of the Fourier coefficients (those corresponding to my frequency band of interest) for computing the covariance matrix, I get a super ill-conditioned matrix.

Does anybody have an idea about how to get around this (that does not involve regularization of the matrix)? or may any other tips on how to approach the problem from a different perspective? Thanks! :)

Well, when you filter the signal to a specific fequency band in the frequency domain, you project your signal to a subspace of $C^T$ . Hence you obtain the a covariance with rank less than $T$.


What is the meaning of inverse signal bandwidth? I have been studying the two ray model in communication systems.

The text I'm studying has two signals, $u(t)$, which is Line of Sight and $u(t-t')$, which is the reflected signal.

It quotes,

  If the transmitted signal is narrowband relative to the delay spread ($B \ll \frac{1}{t'}$) then $u(t)$ is approximately equal to $u(t-t')$.

What is the meaning of this statement?
How is the delay related to the bandwidth?

Sorry if the question wasn't clear the first time.
For the source, one can refer to https://en.wikipedia.org/wiki/Two-ray_ground-reflection_model or Andrea Goldsmith's Wireless Communication (mine is page 34, chapter 2) which quotes the above.

Inverse signal bandwidth is ($\frac{1}{B}$) and inverse delay spread is ($\frac{1}{t'}$). So either ($B \ll \frac{1}{t'}$) or equivalently ($t' \ll \frac{1}{B}$).

So let me rephrase. How can I show that $u(t)$ is approximately equal to $u(t-t')$ given either ($B \ll \frac{1}{t'}$) or equivalently ($t' \ll \frac{1}{B}$).  In other words, what is the proof of the sentence in quotes (or highlights)?

OK, now it's more clear. In a two ray model, you have tow versions of the signal; the line of sight $u(t)$, and the reflected version $u(t-\tau)$, where $\tau$ is the delay spread of the channel.

If $\tau<<\frac{1}{B_s} \simeq T_s$, where $B_s$ and $T_s$ are the bandwidth and duration of the signal, then the shift in the reflected version would be very small compared to the symbol duration. So, the LOS version would be received at time $t_0$ up to $t_0+T_s$, and the reflected version would be received at time $t_0+\tau$ up to $t_0+T_s+\tau$. But $\tau<<T_s$, and thus the two versions are received very close to each others. This answers your first question.

As for the second question, the maximum delay spread of the channel is related to something called the coherence bandwidth of the channel. Roughly speaking, the coherence bandwidth $B_c$ is inversely proportional to the maximum delay spread, in this case $\tau$. If the signal bandwidth $B_s\simeq \frac{1}{T_s}<B_c\simeq \frac{1}{\tau}$, then the channel is said to be frequency-flat fading, which means that the channel frequency response remains contact withing the signal bandwidth. In the time domain, this is translated to interference-free channel, since different symbols don't interfere with each others. On the other hand, if $B_s > B_c$, then the channel is called frequency-selective channel, as different frequency components of the signal are affected differently by the channel. In the time domain, this causes inter-symbol interference (ISI), since the maximum delay spread is larger that the symbol duration, and different symbols interfere with each others at the receiver.


Frequency Domain Filtering Can you create a zero phase IIR filter by transforming its impulse response into frequency domain and only taking the magnitude of that frequency repsone?

I have seen this in an open source project in the following steps:

get filter response
transform it into frequency domain
take the magnitude of the frequency data (-> phase response becomes 0)
multiply with frequency data of a (to be filtered) signal

Does this make any sense? Is this method applicable? I have never seen this described in any literature.

Thank you in advance

This is just "faking" the magnitude response of an IIR filter. The output's magnitude spectrum looks just like it has been filtered by the IIR filter with the given frequency response. Although it may somehow work, there are some limitations:

Frequency-domain filtering is usually much more computationally demanding.
It is not for real-time.

The problem with IIR filters is that they don't have linear phase and cause phase distortion. However, when the filtering is not in real-time (such as when dealing with a stored signal on a computer) there are some alternative non-causal approaches to make the IIR filter zero-phase. The most notable is forward-backward filtering which Matlab's filtfilt also uses. It actually is simple: to make the phase impact symmetric, pass the signal twice from the IIR filter, in opposite directions. After the signal passed the filter (the forward pass), time reverse it and then pass the result again through the filter. So the second (backward) pass cancels the phase impact of the forward pass. However, the magnitude response of the overall filer is the square of the frequency response of the IIR filter.


Frequency and phase range of the FFT of an image I have some questions:

1. What is the horizontal and vertical frequency range (and the steps) of
    the FFT amplitudes of an image?

2. What is the relevant axes range for the phase image?

3. Which sector is physically relevant?

May be you can explain it for the following image:

of which the FFT amplitude image is:

The FFT calculates the discrete Fourier transform (DFT) of an image.
Since it is a separable transform (calculated on rows and columns independently) the explanation is the same for 1D signals.

The DFT of a signal, $X[k]$, is the projection of the signal $x[n]$ onto a complex frequency: $X[k]=\sum_n x[n]e^{-2\pi jnk/N}$.

$k/N$ can be considered to be the frequency. Since $k=0...N-1$ it starts with the lowest frequency (DC or the mean of the signal, $e^0$) and gets the highest frequency - a basis function $[-1,1,-1,1,-1,1,...]$, on the central frequency, $k=N/2$: $e^{2\pi j \frac{1}{2}n }=e^{j\pi n}$. In your example the result has been shifted (using fftshift) so that the lowest frequency is in the middle of the image, and the highest frequencies are in the border.

So, the range is the spatial frequency $0..N-1$ for each dimension. Since the result is complex you can present it using a magnitude image and a phase image, but the axes are the same. A certain point in the 2D transform image $(u,v)$ represents a frequency $e^{2\pi j(xu/N+vy/M)}=e^{2\pi j ux/N}e^{2\pi j vy/M}$ which turns out to be a spatial frequency with a certain direction and "speed".

E.g., $X[u,v]$ where $u=5, v=10, N=M=80$ represents the projection of the image onto the following basis function (here shown is the real part):

In your image, you can see most of the frequencies are concentrated near the origin (i.e. low frequencies), because most of the image is flat. Since there are edges (high frequencies) you can see some high magnitude lines on the transform image.
You can compare it to an image with high frequency content, such as white noise, in which there will be much more high frequency content.


Applying Lowpass (LPF) and Highpass (HPF) Filters to an Image in Frequency Domain in MATLAB What I want is simply to apply low pass and highpass filters to an image in the frequency domain. Does MATLAB image processing toolbox have any commands for this?

You can use spatial Methods using Convolution (conv2 and imfilter).
You can do it in Frequency Domain using the Convolution Theorem property using fft.

Usually, In Image Processing, all filters are built using Gaussian Kernel.
High Pass is basically subtracting Low Pass filtered image from the original image.

You can create the Gaussian Kernel manually or using fspecial.

For fast implementation, have a look here - Fast Gaussian Blur.

P.S.
For people who use conv, be aware that conv is using conv2 and you can reduce the overhead by using it directly.


problem with decimate in scipy 0.18 I use scipy.signal.decimate  to downsample a large spectroscopic data-set.
This has worked superbly so far.

yy = decimate(data.row(i).buffer, n2, ftype = "fir")

I used successively to downsample the spectrum to several lower resolutions.

When I moved from scipy 0.17 to 0.18 (actually 0.18.1) during of a python3 port, a problem appeared.
The first downsampling works, on all the rows; but the next one crashes with the following message:

File "..../lib/python3.5/site-packages/scipy/signal/signaltools.py", line 3049, in decimate
raise TypeError("q must be an integer")
TypeError: q must be an integer

I looked to the code, but with no avail. I also tried the new zero_phase=True argument, but did not help. Needless to say that my n2 argument above IS an integer, the second time as well as the first time!

Any pointers ?

I have found the bug !
n2 is indeed a integer in my case, but its not an int()

print(type(n2)) gives:   <class 'numpy.int64'>   !!!

In the previous version of scipy, it was already a numpy.int64 but it was working correctly.

Maybe a test slightly too stringent...

Thanks for helping everyone.


When I ask Matlab for a zero-padded FFT, does it (or FFTW) ever explicitly create a zero-padded input vector? If I have a length $N$ signal and ask Matlab for an $M$-point FFT, $M \gg N$, of course the computer has to allocate memory for $M$ output samples, but my question is if it ever has to allocate memory for $M$ input and fill it with the input vector along with zeros? Or are FFT algorithms capable of generating the $M$-point result without having to explicitly store the $M$-length padded input?

(Motivation I have a library that can automatically give me zero-padded data (or unpadded data), which I then perform FFTs on. Assume everything before the result of this library call has already happened, and we only care about what we do once those results (padded or unpadded) are available. If Matlab or FFTW don't need to explicitly allocate $M$-long zero-padded inputs to the FFT, then I'd be wasting memory by asking the library to zero-pad for me, and should just ask it for unpadded data.)

There is an FFT trick to perform a $2N$-point FFT on $2N$ real-only input samples using a single complex $N$-point FFT operation. But I've never heard of performing an $M$-point FFT on $N$ input samples without zero padding the input sequence out to length $M$. Ahmed, if you can develop a way to do that you would become famous and someone would write a song about you.


What are F0 counters in Speech signals? What is represented by F0? I can't find any dedicated definition for Frequency F0 counters. Can anyone tell me what does F0 exactly mean, I am confused about that.

Check this source

  To begin with, let’s remember what the fundamental frequency is and in
  what tasks it may be needed. The fundamental frequency, which is also
  referred to as F0, is the vibration frequency of the ligaments when
  pronouncing voiced sounds. When pronouncing unvoiced sounds, for
  example, by whispering or uttering hissing and whistling sounds, the
  ligaments do not vibrate, which means that this characteristic is not
  relevant.
  *Note that the split on voiced/unvoiced sounds is not equivalent to the split on vowels and consonants.

  The variability of the fundamental frequency is quite large, and it
  can vary greatly not only between people (for lower average male
  voices the frequency is 70–200 Hz, and for women it can reach 400 Hz),
  but also for one person, especially in emotional speech.

  F0 is used in a wide range of solutions:

Emotion Recognition
Sex determination (male/female voices)
Speaker deterioration, or splitting the speech into phrases
In healthcare, detection of the pathological characteristics of the voice (for example, using the acoustic parameters of Jitter and

  Shimmer): F0 could be used for the detection of signs of Parkinson’s
  disease 1; Jitter and Shimmer can also be used for Emotion
  Recognition


Determing the frequency of a signal with a sample at a time How to determine the frequency of a signal with a value of the sample and value of the time instant the sample was taken ?

A single sample tells almost nothing about the frequency of a signal.  The only hypothesis that is eliminated is a DC signal of a different value from that of the sample.

The more samples one accumulates over a longer period of time, the better one might be able to fit a maximum likelihood estimate of higher degree models of what the spectrum of the signal might be to those samples, assuming non-aliased samples of a high enough density relative to the bandwidth the signal.


Calculate impulse response when output contributes in input I have an exercise in which I need to find the impulse response for this given system:
$$y(n)=\frac{1}{2}y(n−1)+x(n−1)+x(n)$$
As per my knowledge, I need to find the homogenous solution. My homogenous solution gives me $$y(n)=\frac{1}{2^n} u(n)$$ But the actual answer shows $$y(n)=\frac{1}{2^n} u(n) + \frac{1}{2^{n-1}}u(n-1)$$ I don’t get where the second part came from.
How to solve this?
Thanks!

It should be pretty clear that just computing the root of the characteristic equation can't be enough to find the impulse response. Because in that case the impulse responses of all systems described by
$$y[n]=y[n-1]+\sum_{k=0}^Nb_kx[n-k]\tag{1}$$
would be the same, irrespective of the values of the coefficients $b_k$.
For the given example, the most straightforward (and most instructive) way to determine the impulse response is simply to go through the given iteration until a pattern emerges:
$$\begin{align}h[0]&=1\\h[1]&=\frac32\\h[2]&=\frac32\cdot\frac12\\h[3]&=\frac32\left(\frac12\right)^2\\\vdots\end{align}$$
Can you see any pattern?


Why does decimation make a system time variant? On Wikipedia I read this : "The Discrete Wavelet Transform, often used in modern signal processing, is time variant because it makes use of the decimation operation."

Why does decimation makes system time variant?

HINT:
Think about the following situation: you have an input signal with alternating ones and zeros, and you decimate by a factor of $2$. Now shift the input signal by one sample and compute a new output signal by decimating by a factor of $2$. Is the output signal a shifted version of the previous output signal? If the system is time-invariant, shifting the input signal causes the output signal to shift by the same amount.


Finding specific spots in an audio signal In theory, how would one process an audio signal to find the spots marked in the linked picture below.

Mel Spectrogram

Unfortunately I am not really conform with signal processing and could not find any fitting solutions while looking for similar problems.
I have read a lot about signal processing but might have looked in the wrong places.

Currently, my idea would be: windowing over the signal, applying FFT to each window and check if the frequency/signal strength exceeds a certain value. I am not quite sure if this would work out, plus it seems rather ineficcient.

Can anyone point me in the direction of where to begin and how to identify these spots?

What you are describing is exactly how to plot a spectrogram. Finding the spots can be done with a classic peak-finding function in any software.

If you are looking for coding the algorithm from scratch, here is a method, which is the most intuitive but not the fastest in term of calculations :

Consider your spectrogram as a $n\times m$ array of numbers, with each number position referenced as $(i,j)$. You have to find which points are greater or equal to all of their neighbors.

In pseudo-code it would give something like :

for i from 1 to n
  for j from 1 to m
     if :
       number(i,j) > number (i-1,j-1)
   and number(i,j) > number (i-1,j)
   and number(i,j) > number (i-1,j+1)
   and number(i,j) > number (i,j-1)
   and number(i,j) > number (i,j+1)
   and number(i,j) > number (i+1,j-1)
   and number(i,j) > number (i+1,j)
   and number(i,j) > number (i+1,j+1)
     then :
       save (i,j) as a peak value
   end if
end for
end for

Of course, this is not very sharp since it does not consider whether you are finding a "big" peak or just a smaller one. In that case, you can improve it by setting a minimum offset value like :

if number(i,j) > number(i-1,j-1) + offset

Note that this is just a simple example of a peak finding algorithm and I didn't took the boundary points specific case into account. More efficient methods exist but I can detail each one here.


Uncertainty principle - Duration bandwidth principle In one of NPTEL courses about time-frequency analysis, the professor said that the duration bandwidth principle is $\sigma_t^2 \sigma_\omega^2 \ge \frac{1}{4}$.

He added that the formula making use of time resolution and frequency resolution is a false one. The time and frequency resolutions correspond  to the time "distance" between two samples of the signal and frequency resolution corresponds to the "distance" between two samples in the frequency domain (two successive samples of the Fourier transform of the time domain signal). Can anyone, please, clarify this to me as I get to see the one using time and frequency resolutions in so many papers?

NB: Here is the link to the course I'm talking about http://nptel.ac.in/courses/103106114/

An important  theorem, known as Weyl's, 1931, is:
if  function $s(t)$ and  related functions $ts(t)$, $s'(t)$ are in $L^2$ (square integrable) with the related $\|\cdot\|$ $L_2$ norm symbol then:
$$ \| s(t) \|^2 \leq 2\| ts(t) \| \| s'(t) \|$$
Equality is attained when  $s(t)$ is a modulated Gaussian/Gabor elementary function defined as:
$$ s'(t) / s(t)  \propto t $$
or practically as:
$$s(t) = C \exp [-\alpha(t - t_m)^2 + \imath 2 \pi \nu_m (t - t_m)]$$
found by integration by part + Cauchy–Bunyakovsky–Schwarz.
If one defines  time or frequency location, as a center of mass related to energy as:
$$ E = \int |s(t)|^2 dt = \int |S(f)|^2 df$$
and
$$\overline{t} = 1/E \int t|s(t)|^2 dt \qquad \overline{f} = 1/E \int f|S(f)|^2 df$$
and  energy dispersion as:
$$\Delta t = \sqrt{1/E \int (t - \overline{t})^2 |s(t)|^2 dt }$$
$$\Delta f = \sqrt{1/E \int (f - \overline{f})^2 |S(f)|^2 df }$$
then for  finite-energy every signal  $s(t)$,  with $\Delta t$  and $\Delta f$ finite, one gets:
$$\Delta t\Delta f \geq \frac{1}{4\pi}$$
The limit is attained for some Gaussian variations. My interpretation of the question is is more about the term 'time/frequency dispersion' than the fuzzy concept of 'time/frequency resolution'.


Finding a change in a signal I am looking to some C algorithm(for MCU), where I have these signals :

First graph you can see that there is an event in the middle and it goes back to normal, second and third there is some event at a certain point.
Any of them can happen per session(= when I turn the device there will be 1 behavior)

I have to find  real time when the event begin and ends, and it has to be FAST.

All maximums in all 3 graphs are more or less the same (with 2% diff)
I don't know the values before hand but I can calibrate the "no event" state.
I need 1 algorithm to all of them.

I can do all sorts of basic stuff with averages and LPF, but I guess there is more clever and generic way to do so (?)

Something that can cover all cases and do what my eyes do - I can just tell where is the event.

You need to be sure about:

Noise power is bounded
The power of the signal step is higher than 4 times that noise power.

If the noise power is bounded and the change in the signal level is higher than 4 times that noise power (from the pictures it seems to be so, I mean, the step is more than twice the max value of noise), why not take the derivative of the signal? something like:

step = abs(x[n] - x[n-1]);
if (step > threshold)
{
  // Event has happened:
  // 1. Take timer value
  event_time = get_time();
  // 2. Update the state machine
  state = (state == START_STATE) ? STOP_STATE : START_STATE;
}

Well, the get_timer() is specific to your MCU and there are better ways to implement a Finite State Machine, but that is the idea. It is fast and fulfills the requirements, I would not go further.


high & low pass correlated dither noise question I read a pretty cool article the other day: Robert Katz - Dither noise Probability Density Explained.
It says that if you have two dice (A and B) that you can roll both dice and then...

Re-roll die A and sum A and B

Re-roll die B and sum A and B

Re-roll die A and sum A and B

repeat to get a low pass filtered triangular noise distribution.

It says that you can modify it for high pass filtered triangle noise by rolling both dice and then...

Re-roll die A and take A - B

Re-roll die B and take B - A

Re-roll die A and take A - B

repeat to get a high pass filtered triangular noise distribution.

What i'm wondering is, what is the right thing to do if you want to do this with more than 2 dice? (going higher order)
For low pass filtered noise with 3+ more dice (which would be more gaussian distributed than triangle), would you only re-roll one die each time, or would you re-roll all BUT one die each time.
I have the same question about the high pass filtered noise with 3+ more dice, but in that case I think i know what to do about the subtraction order...  I think the right thing to do if you have N dice is to sum them all up, but after each "roll" you flip the sign of every die.
Can anyone illuminate the correct way to do this?
Edit: I'm looking to get high and low pass filtering with more dice, but am fine with the distribution changing from triangle to gaussian or other shapes.

To achieve high-pass shaping , the “DC term” must be 0. Note that in the 2-dice high-pass case, every roll appears twice, once with a positive sign and once with a negative sign, which enforces the DC-free condition.  Any scheme involving more than 2 dice would still need to meet this condition.
The triangular PDF condition is satisfied by the fact that each sum or difference combines 2 rectangular PDF’s. I don’t know how you can get a triangular PDF by combining 3 rectangular PDF’s, regardless of how they are weighted.

Bob Adams


Savitzky-Golay filter parameters I am trying to smooth a series of data in order to obtain a continuous function that could represent that given data set.
It came out that the Savitzky-Golay method could be a good way.

Now, I don't know much about smoothing and/or interpolate, but I managed to obtain a good result, at a first glance:

The blue points are my data, and the red line is the result from a SG filter applied with a window = 15 and a polynomial order = 13.
Now, I could try some other combinations of window&order, but more generally I would like to know the principles of the best meeting of these two parameters, especially in my case, where the data have no smooth trend (like a sinusoidal).

Let's call the window length $M$, and $N$ is the order of the polynomial. One important thing to realize is that you use a polynomial of order $N$ to approximate $M$ data points. This means that you use $M$ points to compute $N+1$ polynomial coefficients (an $N^{th}$ order polynomial has $N+1$ coefficients). So if $N+1=M$ then you do not smooth at all but you just use the polynomial to interpolate the data points. In your example this is almost the case, i.e. there is hardly any smoothing. Generally, $N$ is chosen considerably smaller than $M$ to achieve some smoothing (and also for numerical stability). The smaller $N$ compared to $M$, the more smoothing you will get. What remains is the choice of $M$, but as far as I know there are no general guidelines. In practical examples I've seen polynomial orders of $4$ or $5$ with much larger values for $M$. With all this in mind it remains to try several combinations of $M$ and $N$ and choose the one that fits best your application.


A non-lossy aka lossless codec without patent-restrictions for audio? The term non-lossy codec means the coding thing having no compression lost, more here. I need such format to process the audio, without destroying the quality. Wikipedia mentions in the context of lossy compression:

Ogg Vorbis (noted for its lack of patent restrictions)
Speex (noted for its lack of patent restrictions)

but the other wikipedia site for lossless things is open-ended, it mentions some patents and some expirations but not whether some codec has no patent-restrictions. Does there exist a lossless (non-lossy) codec without patent-restrictions or some sort of lists/matrix summing up the patent -things?

Perhaps related

Xiph.Org, the organization behind the multimedia formats and software tools here

FLAC is the de facto standard, particularly, notice from the latter link:

  The specification of the stream format can be implemented by anyone
  without prior permission (Xiph.org reserves the right to set the FLAC
  specification and certify compliance), and that neither the FLAC
  format nor any of the implemented encoding / decoding methods are
  covered by any patent. It also means that the reference implementation
  is free software. The source code for libFLAC and libFLAC++ is
  available under the BSD license, and the sources for flac, metaflac,
  and the plugins are available under the GNU General Public License.

  In
  its stated goals, the FLAC project encourages its developers not to
  implement copy prevention features of any kind.[5]


What do you get applying FFT to Music? I am trying to take a section of a song (approx 1.5 seconds)  44100 sample rate equals sample size of 65536 which is 2**16.

I want to find the frequency of the piano note that is played in the window of time.

I have implemented the four1 function - real data starting at array pos 1 and I data at 2 ...

four1 uses are rewritten fortran function so the array does not start at index 0.

My I data is all 0.0.

I have verified that input data looks good ... small sample:  Idx is index into data array

Idx 22023   R -1876.000000   I 0.000000
Idx 22025   R -1720.000000   I 0.000000
Idx 22027   R -1453.000000   I 0.000000
Idx 22029   R -1203.000000   I 0.000000
Idx 22031   R -1017.000000   I 0.000000
Idx 22033   R -799.000000   I 0.000000
Idx 22035   R -386.000000   I 0.000000
Idx 22037   R 75.000000   I 0.000000
Idx 22039   R 243.000000   I 0.000000
Idx 22041   R 327.000000   I 0.000000
Idx 22043   R 759.000000   I 0.000000
Idx 22045   R 1238.000000   I 0.000000
Idx 22047   R 1317.000000   I 0.000000
Idx 22049   R 1335.000000   I 0.000000
Idx 22051   R 1639.000000   I 0.000000

Upon return -- the array contains the R and I of the freq component.

Looking at the first half (real frequency half) -- I get a range of different R and I values but when I calculate the Magnitude (srq of sum sqrs) --- the magnitudes are almost constant the complete real half of the freq range

Sample of output -- I normalized the Magnitude to the Max Magnitude found

Idx 13841   R -605090688.000000   I -74311704.000000    M 0.958672
Idx 13843   R -605982336.000000   I -74471496.000000    M 0.960094
Idx 13845   R -605715712.000000   I -74230088.000000    M 0.959632
Idx 13847   R -605639744.000000   I -74081136.000000    M 0.959485
Idx 13849   R -605420992.000000   I -74552360.000000    M 0.959234
Idx 13851   R -606319744.000000   I -73444392.000000    M 0.960425
Idx 13853   R -605728512.000000   I -73858336.000000    M 0.959581
Idx 13855   R -605823424.000000   I -73527352.000000    M 0.959666
Idx 13857   R -605311424.000000   I -73650440.000000    M 0.958890
Idx 13859   R -606479360.000000   I -73952432.000000    M 0.960771
Idx 13861   R -606198336.000000   I -74041552.000000    M 0.960349
Idx 13863   R -605836544.000000   I -74178936.000000    M 0.959811
Idx 13865   R -605548864.000000   I -73529936.000000    M 0.959238
Idx 13867   R -605573312.000000   I -73070880.000000    M 0.959189

There is little to no change in Magnitude across freq bins.

Should I not be able to detect the frequency of that piano note played in that time window?

Am I misusing the FFT?

Thanks in advance for any reply.

Yes, you should be able to see the component of the piano note. But there are quite a few pitfalls on the way. I would suggest you experiment a little to get a 'feel' on how fft works first. Use GNU Octave or Matlab to do simple experiments, and then you can plug in your real data, and see if it works. For example, in Octave I can easily generate a sine wave of 440 Hz and do the FFT:

t = linspace(0, 2*pi, 44100);
s = sin(440.0 * t)              % This is the 440 Hz tone
f = abs(fft(s))                 % Take the fft (magnitude)
plot(f)
xlim([0,1000])                  % Show only the 1000 first bins

Then you can read your samples and do the same. At least, this way you
exclude insecurity about the fft function you're using. You could also
save the generated sine, and feed it to your four1 function, and check
if the result is correct.


Determine the time length of audio training samples? Currently I am trying to use MFCC feature for audio training in order to do human voice classification (can later distinguish the person who speaks).

I have problems to determine the proper time length of audio training samples? To be clear, say I have computed 12d-MFCC features over 20ms window (10ms overlapping), if the length of the audio training samples is 1s, this will result in a 12 x 100 feature vector (as there are 100 frames during the 1s time given 10ms shifting time step). In practice, what's the proper time length of audio training samples? Of course, long time length will results in long MFCC feature vector.

It's a very bad idea to stack together consecutive frames in a single vector. By doing so, you are training a classifier to recognize the exact same sound as the one in your training database, pronounced at the exact same speed.

For your application (which appears to be speaker identification), you should consider each individual frame as a training instance (aka "bag-of-frames" models). That is to say, each training sample will give many 12-dimension feature vectors.


simulink discrete integrator block problem I am testing the Simulink discrete integrator block as shown in the picture below. The simulation is fixed step with $T_s = 1e-5$ and the cosine wave block setting and the scope display is also shown below.

For the discrete integrator block, the gain value K is set to be 1 and the integration method was chosen to be "Accumulation: Forward Euler". However, how come the magnitude of the result is so big as shown in the top part of the scope picture? I thought integrating a cosine wave will produce a sine wave with the same magnitude. If the integration method is chosen to be "integrator: Forward Euler", the result becomes very small, as shown below. It just can't produce a sine wave with magnitude of 1. How to get the correct result?

All of a sudden I realized I made a silly mistake.
$$
\int \cos(\omega\cdot t)dt=\frac{1}{\omega}\sin(\omega\cdot t)
$$
which means the "integrator: Forward Euler" choice produced the correct result.


Mitigate floating-point numerical errors for very-low corner frequency low-pass filter with DSP I'm designing a low-pass filter for a digital signal processing application that ideally just passes a very small bandwidth above DC. I'm using an IIR biquad filter for this, where the coefficients are derived using the instructions here. A smaller bandwidth leads to a longer filtering time (larger time constant) but yields a more accurate result whereas a larger bandwidth can be filtered faster but is less accurate. Both of these are valid use cases.
Here's the code I've got
#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz

# calculates filter coefficients using link above
# fc is corner frequency, fs is sample freq
def iir_lp_coeffs(fc, fs):
    w0 = 2 * np.pi * fc / fs
    q = 1 / np.sqrt(2)
    alpha = np.sin(w0) / (2 * q)
    b0 = (1 - np.cos(w0)) / 2
    b1 = 1 - np.cos(w0)
    b2 = b0
    a0 = 1 + alpha
    a1 = -2 * np.cos(w0)
    a2 = 1 - alpha
    b0 /= a0
    b1 /= a0
    b2 /= a0
    a1 /= a0
    a2 /= a0
    a0 /= a0
    return (
        np.array([b0, b1, b2], dtype=np.float64),
        np.array([a0, a1, a2], dtype=np.float64),
    )

fc = 2  # low pass corner frequency (Hz)
fsample = 500e3
b, a = iir_lp_coeffs(fc, fsample)

w, h = freqz(b, a, worN=int(1e6), fs=fsample)
fig, ax = plt.subplots()
ax.plot(w, 20 * np.log10(abs(h)))
ax.set_ylim(-40, 10)
ax.set_xscale("log")
plt.show()
print(w[0:10])
print(abs(h[0:10]))

The current settings use 64-bit floating point with a cutoff frequency of $2\,\text{Hz}$. This all works fine, and I can even decrease the corner frequency substantially as long as I increase the granularity of freqz (with worN=).
For instance here's a plot of the gain response with the above code (note that I've cut the x axis off at the higher frequencies):

However, my actual application requires 32-bit floating point. When I do this (set dtype of iir_lp_coeffs to np.float32), I get non-unity gain in the passband. For instance, here's a gain response with fc=10 using 32-bit:

If I set the corner frequency higher, the gain response looks correct again (e.g., fc=100 looks fine).
Am I running up against the limit of what's possible with 32-bit FP? Or, is there another strategy that would allow me to get away with the lower precision of 32-bit? Have I correctly diagnosed this issue as a floating-point issue?

I think your issue might be coefficient quantization and filter topology. A direct form biquad has poor quantization effects about 0 and π radians. It’s easier to analyze such effects in fixed point, but even though floating point has a much larger range, it still has shortcomings. In particular, if you add a very small number to a very large number, the small number disappears because it can’t be aligned for the operation in the available number of mantissa bits. This can cause the order of operations to affect the result. For instance, where S is a small number and L is large, L - L + S = S, but L + S - L = 0.
Udo Zolzer covers the differences between several filter structures in his book Digital Audio Signal Processing. I borrow direct form quantization effects on pole locations from the book:

See how precision is lost near 0 and π. Other filter topologies might have higher precision near 0, while being much worse near π, which may be a good tradeoff for uses such as yours. The Gold and Rader form has a very even distribution, it looks like a perfect grid.
Another simple and popular filter that has good quantization characteristics at low frequencies is the "Chamberlin" state variable filter. There are improved versions of this filter, as it has problems at higher frequencies (from about one-sixth of the sample rate and up), but the plain Chamberlin is very good at low frequencies, where you need it.
See my article on the Chamberlin state variable filter here:
The digital state variable filter
Zolzer presents modified Chamberlin structures here:
The Modified Chamberlin and Zölzer Filter Structures
In particular, see the graph of quantization effect near zero for the Chamberlin structure—very dense near zero, at the expense of poorer performance at high frequencies, compared to the direct form graph:


Evaluate performance of FFT Windows lately i've been interested in FFT windowing, i've experimented various windows in Matlab, that made me more curious:

There are normalized parameters (CG, Scalloping loss etc...) for every window

How those parameters were measured? for what type of signals? what are the frequency, amplitude of those signal? what is the Frequency sampling used? what is the length N of the  FFT?

is there a standardized framework to use for evaluating windows.? what if i want to do an evaluation my self, what initial conditions i should choose ?

The good thing about filters is that they are linear. Therefore you do not need to specify the signal when you are analyzing the filter. You just analyze its impulse response (i.e., filtering an ideal Dirac pulse) and the Fourier transform of that (the transfer function). Based on this analysis, you can judge what the filter would do to any given signal.

For windowing, almost the same is true, just a little backwards: windowing means multiplication in time domain, which is convolution in frequency domain. So it's almost the same as filtering, just exchanging the role of time and frequency domain. Therefore, we can do the same: analyze the window function and its Fourier transform itself, not assuming any specific signal. That's what would happen for a single Dirac pulse in frequency (i.e., a constant DC signal). Then, if you want to know what happens when you window a given signal, just apply convolution in frequency domain (e.g., for a pure tone, that's just shifting it).

That said, most of the table can be filled by computing the Fourier transform of the window function and then taking the corresponding measurements: how large/wide are the mainlobes, etc. If you cannot compute the Fourier transform analytically, you can try to approximate it, using a Discrete Fourier transform. In this case, use a large $N$, to make sure your estimation errors remain negligible.

*edit: Regarding your question on the estimation of power: if you want peak heights to be consistent you'll have to normalize your windows accordingly. Since the spectrum of a windowed signal is $X_h(f) = X(f) * H(f)$, i.e., convolution with the spectrum of the window function, the normalization of $H(f)$ matters. For instance, for a small peak to come out in the correct height, you might want to have $H(0)=0$, which means that the coefficients in time domain should sum to one. Alternatively, if you want your windowed signal to behave like the unwindowed one, take into account that applying "no window" means implicitly applying a rectangular window of amplitude 1, so that its coefficients sum to $N$.

In Matlab:

clear;
N = 999;
fs = 2000;
f0 = 200;

taxis = (0:N-1)/fs;
faxis = (0:N-1)/N*fs;
sig = sin(2*pi*taxis*f0);

win = flattopwin(N).';
sig_win = sig.*win;

win_normalized = win/sum(win)*N;
sig_win_normalized = sig.*win_normalized;

figure(1);
clf;
plot(faxis,abs(fft(sig))/N);
hold on;
plot(faxis,abs(fft(sig_win))/N);
plot(faxis,abs(fft(sig_win_normalized))/N);
legend('Rect','Flattop','Flattop normalized')
xlabel('f [Hz]');

set(gca,'YScale','log');

Result (zoomed in around 200):

You can see how with the normalization you can make the flattop window behave like the rectangular window.

Btw: the units of this are not dBm. It's a discrete Fourier transform, transforming Volts into Volts (or unitless quantities into unitless). If you want dBm you have to properly normalize to make sure your FFT approximates a Fourier transform (it involves the sampling time $t_0$) and you'll have to estimate power instead of amplitude.


Confusion in Phase Shift Keying What is the difference between Digital phase modulation(DPM) and Phase shift keying (PSK)? http://suraj.lums.edu.pk/zartash/publications/2007-DigitalPhaseMod.pdf It says that PSK is a special case of DP.

Phase shift keying is linear modulation.

Digital phase modulation is not necessarily linear.